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SUMMARY 


This study was composed of three parts. The first part 
involved the extension of an existing transient two-dimen- 
sional numerical code for an electrothermal deicer so that 
it would simulate the situation where a variable thickness 
ice layer existed at the outer surface. The Enthalpy Method 
was used to simulate the phase change, and Gauss-Seidel it- 
eration was used to solve the resulting system of finite 
difference equations. A set of criteria were developed for 
determining when a variable thickness ice layer had an ef- 
fect on deicer performance. 

The second part of this study was the acquisition and 
analysis of experimental data. The test model consisted of 
a section of a Bell UH-1H helicopter blade equipped with an 
electrothermal deicer manufactured by the B. F. Goodrich 
Company. A total of fifty-two thermocouples were utilized 
to document the thermal response of the blade and deicer as- 
sembly. The tests were conducted in the Icing Research Tun- 
nel at the NASA Lewis Research Center, and consisted of four 
phases; dry air tests; wet air tests; ice accretion tests; 
and deicing tests. A total of two hundred and eleven tests 
were run, from which ten readings, five dry runs and five 
deicing runs, were selected for further analysis. This re- 
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duced set of data were examined for thermocouple response 
consistency and, within experimental error, the thermocouple 
readings were found to be independent of the three cross- 
sections on the blade where measurements were taken. The 
dry run temperature responses at the abrasion shield showed 
test independence when correlated as AT/(AT)avg versus posi- 
tion. This was not true for the deicing runs. For both the 
dry and the deicing runs, the heater response, when corre- 
lated in the same manner, showed test independence for all 
ten runs. This behavior was a physical characteristic of 
the blade. In contrast, the response of the D-spar thermo- 
couples was found to be almost entirely independent of posi- 
tion within each test. In the deicing runs, the experimen- 
tal temperature response data clearly showed when melting, 
shedding or refreezing occurred. These tests illustrated 
that the criterion for shedding in the three cases where it 
did occur was that the abrasion shield interface temperature 
was 32-34°F. 

The third part of this study concerned the validation of 
a one-dimensional transient thermal model of an electrother- 
mal deicer by comparison of the predictions with the experi- 
mental data. The physical properties and the geometry of 
the test model were determined and used as input m the nu- 
merical simulation. Flat plate and cylindrical correlations 
were used to calculate the outer surface convection coeffi- 
cients for both the dry and the deicing runs. The Enthalpy 
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Method was found to effectively model the phase change which 
occurred, and the ice shedding algorithm employed in the 
simulation was also evaluated. In general, the one-dimen- 
sional code showed good comparison with the experimental 
data, with the comparisons being better at the higher free- 
stream temperatures. The code definitely contains suffi- 
cient physical information so as to adequately model the 
thermal response of an electrothermal deicer assembly at po- 
sitions where two-dimensional effects are small. 
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A NUMERICAL AND EXPERIMENTAL INVESTIGATION 


OF 

ELECTROTHERMAL AIRCRAFT DEICING 


PART I 

INTRODUCTION 

Ice accretion on aircraft components during flight has 
long been recognized as a serious aviation hazard. The 
presence of the ice decreases lift and increases drag, as 
well as increasing fuel and power requirements. Military 
and commercial aircraft that are not protected cannot fly in 
icing conditions. As a result, there is a continuing need 
for the development of systems that either prevent the ice 
from forming or remove it after accretion. 

Preventing the formation of ice is known as anti-icing, 
while allowing it to form and then removing it is known as 
deicing. There are several ways of achieving either goal. 
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but some methods are better than others depending upon the 
particular application. The main thrust of the present work 
concerns the analysis of electrothermal deicing systems as 
applied to helicopter blades. As will soon become apparent, 
there are some limitations peculiar to this application that 
predicate a certain type of system. But first, a brief 
overview of the different anti-icing and deicing methods 
available is in order. 

There are two commonly used types of anti-icing systems, 
chemical and thermal. Chemical systems use a freezing point 
depressant much like an antifreeze. The chemical is spread 
over the area to be protected, thereby lowering the freezing 
point and preventing ice accretion. This method is limited 
because of the difficulty encountered in maintaining a uni- 
form distribution of the chemical over the surface. There 
is also a reservoir that frequently needs to be refilled, 
and holes used to transmit the antifreeze to the surface may 
become plugged. As a result, chemical systems are used most- 
ly for windshield protection. 

Thermal anti-icing prevents ice accretion by raising the 
surface temperature above the freezing temperature. Gener- 
ally this is accomplished by passing hot compressed air 
through passages beneath the area to be protected. Energy 
requirements are quite high to achieve this. However, jet 
aircraft usually utilize this method because an ample supply 
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of hot compressed air may be bled from the first com- 
pressor stage of the engine. 

Three types of deicing systems worthy of mention are: 
pneumatic boots, electromechanical impulse and electrother- 
mal heating. A pneumatic boot is made of a flexible rubber- 
like material that is laminated to the area to be deiced. 
It is inflated like a ballon with compressed air, which me- 
chanically breaks the ice from the surface. Boots are rela- 
tively simple and efficient, but they increase drag and re- 
quire frequent maintenance to insure reliability. 

Electromechanical impulse is a new method still in devel- 
opment. A series of electromagnets are pulsed m cycles, 
flexing a metal abrasion shield and thereby mechanically 
cracking off the ice. This method appears to be extremely 
energy efficient, but it has not yet been applied to heli- 
copter blades because it is still in the developmental 
stage . 

Electrothermal deicing uses electric heaters laminated to 
the area to be p'rotected. Dynamic forces remove the ice af- 
ter its adhesion to the blade surface is destroyed by melt- 
ing a thin layer of ice next to the surface. The heaters 
may be cycled to reduce energy requirements. This method 
was first applied to deice propellers and lends itself par- 
ticularly well to rotary wing aircraft. 
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Helicopters use electrothermal deicing exclusively to 
protect the blades. Using pneumatic boots on helicopters is 
still in the experimental stage. However, some success has 
been achieved, and boots may become commercially available 
in the near future. Electromechanical impulse is still in 
development and has not been applied to helicopters. At 
this time, electrothermal is the most advantageous method 
available to deice helicopter blades. 

Electrothermal deicers are assembled in a composite con- 
struction. A wide variety of materials may be used but some 
are better suited than others. Baliga [1] and Stallabrass 
[2] have examined the effects of different materials and 
thicknesses. A typical example is shown in Figure 1 and 
consists of: 

1. Substrate - This is actually the aircraft blade and 
can be made of many different types of materials. An 
aluminum alloy is considered here; 

2. An inner layer of insulation - This is necessary to 
provide both thermal and electrical insulation. Re- 
sin impregnated cloth is commonly used for this pur- 
pose. A thickness of 0.01 to 0.02 inches is neces- 
sary to fulfill electrical insulation requirements; 

3. Heater element - Woven mats or metal ribbons are 
used. Woven mats have thicknesses as great as 0.02 
inches. Ribbons range from 0.001 to 0.005 inches m 
thickness. The individual elements are 0.5 to 1.0 
inches wide; 
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4. Outer insulation 


Good electrical insulators are 


usually poor conductors of heat. As a result the in- 
ner insulation has to be thicker than the outer insu- 
lation in order to direct most of the heat outward. 
A ratio of at least 2:1 has been recommended by Stal- 
labrass; and 

5. Abrasion Shield - This is necessary to protect the 
heater and insulation from the environment. Stain- 
less steel is often used, ranging from 0.01 to 0.02 
inches in thickness. Its high thermal conductivi ty 
serves to transmit heat laterally and reduce the ef- 
fect of nonuniformities m heater output. 

A major difficulty with this type of construction results 
from imperfect bonding between layers. Great care must be 
taken during the lamination process in order to prevent air 
gaps. However, small gaps usually do exist which inhibit 
the rate at which heat is transferred and, therefore, impede 
deicer performance. 

Another difficulty that decreases performance is nonuni- 
formity in heater output. Gaps exist between the heater el- 
ements in order to provide electrical insulation. These 
gaps are about 0.08 inches wide for woven mats and 0.04 
inches wide for metal ribbons. These heater gaps cause hot 
and cold spots to occur at the abrasion shield-ice inter- 
face, which delays melting and subsequent shedding of the 
ice . 
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There are many factors that affect the design of elec- 
trothermal deicers. A combination of experimental data and 
analytical analysis is thought to be necessary in order to 
understand the physics occurring during the deicing process, 
and thus to be able to make recomendations concerning deicer 
design . 

The first objective of the present study involved the in- 
corporation of a variable ice thickness algorithm into an 
existing computer model. The second objective was the ob- 
tainment of experimental deicing data using the NASA Lewis 
Icing Research Tunnel. Finally, the third objective was to 
compare the experimental data with the heat transfer models 
developed by Baliga [1], Marano [3], and Chao [4] in order 
to determine the value of these models for use in electroth- 
ermal deicer design. 
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PART II 


A SOLUTION TO THE VARIABLE ICE THICKNESS PROBLEM 

Chapter I 
LITERATURE REVIEW 

1.1 ANALYTICAL HISTORY 

An analytical approach to the phase change problem was 
first proposed by Stefan (1889). Stefan dealt with the 
transition from solid to liquid or liquid to solid. These 
types of phase change problems are classically referred to 
as Stefan problems but are also known as ablation problems. 

Arpaci [5] has formulated the classical ablation problem 
which involves a nonlinear boundary condition at the moving 
solid-liquid interface. This nonlinear equation at the in- 
terface complicates the problem immensely and makes an ana- 
lytical solution very difficult. Reference [6] presents an 
extensive review of most of the analytical and numerical 
techniques used for phase change analysis. 

Predicting temperatures in composite bodies is almost as 
formidable as the ablation problem. A Laplace Transform 
technique proposed by Carslaw and Jaeger [7] is effective 
for one or two layers. Finding the inverse transform is 
difficult for more than two layers. Campbell [8] proposed a 
method that uses an analogy between heat flow and the trans- 
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mission of electricity along a power line. Stallabrass [2] 
incorporated this method m verifying the accuracy of his 
numerical technique. 

Chao [4] provides a more complete summary of the analyt- 
ical techniques available for ablation problems and compos- 
ite bodies. All of these methods are rendered impractical 
by the large number of calculations required for each temp- 
erature. The capacity of computers for executing large num- 
bers of calculations inexpensively and efficiently makes 
them ideal tools for solving phase change and composite body 
problems . 

1.2 NUMERICAL HISTORY 

Recent interest in the electrothermal deicer problem has 
generated several numerical models. One-dimensional models 
were proposed by Stallabrass [2], Baliga [1] and Marano [3], 
and two-dimensional models have been developed by Stalla- 
brass [2] and Chao [4]. All of the models consider phase 
change and use finite difference techniques. 

Finite differencing is an approximate technique for solv- 
ing boundary value problems. Carnahan, Luther and Wilkes 
[9] discuss the various differencing methods that are com- 
monly used and their advantages and disadvantages. The fi- 
nite difference method discretizes the continuous time and 
space domains into a grid of nodes. A system of difference 
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equations based on this grid replaces the governing differ- 
ential equations and boundary conditions. This system of 
equations can be solved for a particular node at any point 
in time and space. Accuracy depends upon the particular so- 
lution method used and the grid spacing. 

There are two basic classes of finite difference tech- 
niques: explicit and implicit. The explicit method uses 
only nodal values from the previous time step in the finite 
difference equations. A distinct advantage to this method 
is that the resulting set of equations can be solved direct- 
ly for the nodal values at the present time step, without 
the use of inversion or iteration procedures. This reduces 
the total number of computations and saves on computer time. 
Explicit methods do not, however, guarantee convergence. 

The increments in time and space must satisfy the following 

2 

equation: aAt/(Ax) <1/2, were a is the thermal diffusivity, 
and At and Ax are the time and space increments, respective- 
ly. This convergence criterion forces a time step of 0.001 
seconds or less, which increases the number of calculations 
needed. Stallabrass [2] and Gent and Cansdale [10] used an 
explicit algorithm scheme which was forward differenced m 
time and central differenced in space. The truncation error 
for this method was first order in time and second order in 
space . 
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Implicit methods use nodal values from the present time 
step as well as from the previous time step. This forms a 
system of equations that must be either inverted or solved 
by iteration at each time step. An advantage of the implic- 
it method is its unconditional stability, which allows a 
time step as large as 0.1 seconds. The governing energy 
equation for electrothermal analysis is a parabolic partial 
differential equation. Von Rosenberg [11] indicated that 
the Crank-Nicolson implicit method is the preferred method 
for this type of differential equation. The truncation er- 
ror for this method is second order correct in both time and 
space. Baliga [1], Marano [3] and Chao [4] all employed the 
Crank-Nicolson implicit finite difference technique m their 
simulations . 

Several numerical methods have been proposed for handling 
the phase change. Baliga [1] used a technique suggested by 
Bonacina et al. (12). This approach modeled the latent heat 
effect with a large change in heat capacity over a small 
temperature interval around the melting point. The thermal 
conductivity was varied linearly over the interval. Stalla- 
brass [2] accounted for the phase change by holding a node 
at the melting point until enough energy had been trans- 
ferred to completely melt the nodal volume. Marano [3] and 
Chao [4] used the enthalpy method which is also known as the 
weak solution method. These methods to simulate the phase 
change are comparable, but the enthalpy method is preferred 
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because of its formalism and the eaae with which it is ap- 
plied. The enthalpy method is used in the present study and 
has been employed m a variety of problems in the recent 
literature [6,13,14,15,16]. 

The enthalpy method uses the conservation of energy equa- 
tion formulated m terms of two dependent variables, temper- 
ature and enthalpy. Predicting the location of the solid- 
liquid interface is not required because this is determined 
by nodal enthalpy alone. The temperature at any node can be 
calculated from the known enthalpy- temperature relationship. 
This method requires the use of Gauss-Seidel iteration to 
solve the resulting set of difference equations. The equiv- 
alence of this method and the classical formulation of the 
ablation problem was proven by Atthey ]14]. 

Temperatures predicted by the enthalpy method tend to os- 
cillate about the true values once melting has started. 
This occurs because a node is held at the melting point for 
a finite period of time until the nodal volume melts. The 
techniques used by Stallabrass [2] and Baliga [1] also ex- 
hibit this phenomenon. Voller and Cross [13,16] have point- 
ed out that this leads to unrealistic results and they have 
developed a criterion for deriving the true temperatures 
from the oscillating values using the nodal enthalpy. This 
criterion is discussed in the Appendix. 
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Chapter 2 

FORMULATION OF THE NUMERICAL METHOD 


2.1 METHOD OF APPROXIMATING THE VARIABLE ICE THICKNESS 

The previously mentioned two-dimensional numerical simu- 
lation code developed by Chao [4] has the restriction of a 
constant ice layer thickness. In order to simulate blade 
sections that have variable ice layer thicknesses, a stair- 
case approach was used in the present study. 

This staircase approach approximates the ice shape with a 
series of discrete rectangular nodes. The x and y direc- 
tions are discretized into dx and dy segments. If the lce- 
ambient interface passes above the centroid of a node, the 
node is added to the ice shape. If the interface passes be- 
low the centroid of a node, the node is not included. In 
this manner', a rough approximation of the ice shape is 
formed, as shown m Figure 2. This approximation may be im- 
proved by further decreasing the sizes of dx and dy. Ice 
accretion shapes usually have a rough surface that can be 
approximated rather naturally using this approach. It is 
also simple and easy to apply. 
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2.2 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


The assumptions made m formulating the mathematical mod- 
el are: 

1. The thermal and physical properties of the material 
composing each layer of the blade are different, but 
do not depend on temperature ■ 

2. The ambient temperature, blade air temperature and 
all heat transfer coefficients are constant; 

3. There is perfect thermal contact between each layer; 
and 

4. The density change due to melting is negligible. In 
other words, the effect of the volume contraction of 
the ice as it melts is neglected. 


Using these assumptions, the governing two-dimensional 
heat transfer equation can be reduced to: 


3T 


P J C pJ 3t 


3 2 T 


3 2 T, 


3y 


2 > + ‘‘j 


(1) 


where j represents the layer in question, as given by 
3 = 1: blade or substrate 
3 = 2: lower or inner insulation 
3=3: heater 

3=4: upper or outer insulation 
j=5: abrasion shield 
and where 

p, = density of the j^layer; 
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A 4-U 

C . = specific heat capacity of the j layer; 

P3 

Tj = temperature in the 3 ^ layer; 

Kj = thermal conductivity of the 3 layer; 

q^ = rate of heat generation per unit volume of the 3 

layer; 

t = time variable; and 
x,y = space coordinates. 


th 


The heat generation term, q_^ , is zero for all layers except 
for the heater, where it becomes a function of time, q^ (t) . 


A slightly different approach must be used for the ice 
layer (3 = 6). The physical properties change as the ice 
melts. Because of this, the thermal conductivity, K^, must 
be included within the first partial time derivative. The 
governing differential energy equation becomes: 



( 2 ) 


At this point, it is advantageous to replace p c C C T, with 

D pD D 

H , , representing the enthalpy at positions within the ice 

D 

layer. Carrying out this operation on equation ( 2 ): 


3H, a 3T t a 3T a 

_A =. i — (v _A . L_ 6 

at 3x 3x > ay '*6 3y 


where 
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H = enthalpy per unit volume within the ice-water 
6 

layer; 

T = temperature within the ice-water layer; and 
6 

K = thermal conductivity within the ice-water layer. 

6 

This equation is the governing equation for the enthalpy 
method. It can be used to solve for the temperature at any 
point within the ice layer. Enthalpy can be related to 
temperature using the following relationships: 


H 


6 


C 

ps 


< T 

D 


P 1 C pl <T 6 ' V 


+ 0. 


(C T 
ps H 


* V* T 6 


> T 


(4) 


where 


D , C = physical properties of the solid phase; 
s ps 

A 

, C ^ = physical properties of the liquid phase; 


m 


= melting temperature; and 


= latent heat of fusion per unit mass. 
Inverting equation (4) yields: 


H /p c 

o s ps 


, H. < H 
b — so 


, H < H. < H. 
sm 6 lm 


( V H i» )/p i c P i + t »' h 6 i h i d 


(5) 


with 


In 


PC T 
s ps n 


P. (C T + L.) 
l ps m f 


( 6 ) 
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where H gm and H are the melting point enthalpies of ice 
and water, respectively. 


The boundary conditions at the layer interfaces, inner 
ambient surface, and outer ambient surface are listed below. 
1. At the interfaces between the layers, the temperature 
and heat fluxes are continuous: 


Mli 

3T 


j+l|l 


-K 


J 3 y 


- -K 


j+l 3y 


+_1 

1+1 


j - 1 5 


J ” 1 5 


(7) 


where I denotes the interface. 

2. Newton's law of cooling is used at the inner surface- 
ambient and outer surface-ambient interfaces. For 
the lower or inner boundary: 


3T 


J 


V <T j 1 


T bl ) • J " 1 


( 8 ) 


where 

h^ = lower boundary convective heat transfer coeffi- 
cient ; 

T, , = lower ambient (blade air) temperature; 
bi 

and the subscript, 1, denotes the lower surface. 
Correspondingly, for the upper or outer boundary: 

2 ’ h b2 (T j| 2 - T b2>* J ‘ 6 

where 

h = upper boundary heat transfer coefficient; 
b2 



16 



T ^2 = upper surface (external flow) ambient 

temperature; 

and the subscript 2 denotes the upper surface. 


If there is a variable thickness ice shape at the 
upper surface, the boundary conditions are much more 
involved. There are actually three different bound- 
ary conditions possible depending upon the shape of 
the ice-ambient interface. 

For a horizontal surface: 



ay 


2 


h b* < t 6 2 " T b2) 


(9a) 


For a right-facing vertical surface: 


3 t 6 
- k* — 
ax 


Z 


h b2 < t 6| 2 " T b2) 


For a left-facing vertical surface: 


(9b) 


3*6 

k 6 

3x 


z 


h b2 < t 6| ^ ” T b2) 


( 9c ) 


The constants in equation (9) are defined as 

h = upper surface convective heat transfer coeffi- 
b2 

cient ; 

T = upper air (external flow) ambient temperature; 
b2 
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and the subscript, 2, denotes the upper ambient in- 
terface . 

3. Insulated boundary conditions can be used to repre- 
sent the symmetry at the centerline of the heater el- 
ement and heater gap. 


lil 

.'X 


g.J 


ST 

3x 


h.J 


. J = i. .. , 6 

. J - 1 6 


(10a) 


(10b) 


where the subscript, g, denotes the position of the 
gap centerline and h represents the position of the 
heater centerline. 


The initial temperature distribution can be constant or a 
function of position. In this study, only a constant ini- 
tial temperature distribution was considered. 


These equations must be expressed m terms of dimension- 
less parameters. Those chosen were: 


e = 


ref 


- _ X - _ y 
X = L ’ y = L ’ “j 


°J C PJ 


where 


T _ = the reference temperature (taken to be 32 °F m 
ref 

this study) ; 

L = the total thickness of the composite slab; and 

a = the thermal diffusivity of the j*"* 1 layer. 

j 
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These definitions were substituted into equations (1) 
through (10), resulting in the following set: 

For each layer of the composite blade: 




p c 

3 


T 

PJ 


ref 


(ii) 


For the ice layer: 



oc. 

L 3x ' 


3x 


30, 


(K, 


3y 


3y 


-)] 


(12) 


P C 0, T 
s ps 6 ref 


’ 6 6 < T m 


p, C , T , (0, - 0 ) + p, (C 0 T + 1 ), O > T 
1 pi ref 6 m 1 ps m ref f 6 m 


H,/p C T , 
6 s ps ref 


, 11, < H 
6 ~ sm 


, H < H, < H, 
sm 6 lm 


H,/p. C 1 T ( c 0 + + e . H 1 H 

6 1 pi ref ps m T__ , m b — lm 


Pi 

At the layer interface points: 


ref 


(13) 


(14) 
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J 
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D J+i i 


. j 
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30 
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. J 
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(15a) 


(15b) 


At the lower and upper boundaries: 


30, 


J *y. 


h bi L V b bi ) ’ J = 1 


DO 


-K 


J 3y, 


= h b2 L (n j 2 - °b2 ) • J “ 6 


(16a) 

(16b) 
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Lastly, the insulated boundary conditions at the heater and 

gap centerlines become: 

ao I 


3x 


g.J 


j - i. 


6 


(17a) 


dO 


|h,J 


1 6 


(17b) 


2.3 NUMERICAL TECHNIQUE AND FINITE DIFFERENCE FORMULATION 


Finite difference formulations require that the continu- 
ous time and space domains be discretized into a set of 
nodes on a grid. For this study, there are three dimensions 
to be discretized. The spatial coordinates, x and y, are 
divided into segments, A x and AY- The time coordinate, t, 
is discretized into At. The set of governing equations, 
initial conditions and boundary conditions must be cast into 
finite difference form. This set of equations is solved .for 
each point of the spatial grid at each time step. A central 
difference scheme was employed, resulting m the following 
equations for the first and second partial space deriva- 
tives : 


dE 




a x 


A+LA _ — ichJs + o (ai ) 2 

2 Ax 


a 2 E 


i.k 


a x 2 


^1+1 k ~ 2 ^i k + ^i-1 k -2 

- - ■ » + 0 <Ax) 

(fix) 


where 


(18) 

(19) 


the i and k denote the nodal position on the two-dimension- 
al grid of nodes, and the E represents any desired nodal 
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quantity. These are second-order accurate finite difference 
analogs. The first partial time derivative may be written 
in this form as: 


at 

7T 


iA 


0 

1 




+ o (at) 2 


( 20 ) 


where the superscript A indicates the value from the previ- 
ous time step and the superscript o denotes the value at the 
half time step between the previous and current time levels. 
Equation (20) is commonly known as the Crank-Nicolson formu- 
lation for the time derivative. The most advantageous 
aspect of this formulation is the higher accuracy resulting 
from a second-order truncation error. The spatial Crank-Ni- 
colson analogy can be derived by taking the average of the 
present and previous time levels. The first and second spa- 
tial partial derivatives become: 


,)U 


l.tc 


3E 


iU <- 


iA 


3x 


3x 


3E 

+ — 
j 3* 


iA 


1+1, k 


- E 


+ E, 


* E, 


-.2 


JAA_ — iilA L~ 1 A + n (ax) 

A Ax . 


( 21 ) 


A 

3x 2 


v.k 


3 2 E 


1/2 (- 


iA 


3x 


|j 


3 2 E 

I? 


iA 


v -2E + E + e ^ -2!'^ + f ^ 

b M,k IA 1+1. k IJLJL— kOJj 4 o < A x) 2 

2 (Ax) 2 


( 22 ) 
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Equations (20) through (22) can be substituted into equa- 
tions (11) through (17) to form a set of nodal equations 
that can be solved for any point within the discretized time 
and space grid. 

2.3.1 Finite Difference Equations for the Composite Body 

This study is based on the two-dimensional model devel- 
oped by Chao [4] to study the effect of a gap m the heating 
element. A complete explanation of the derivation of the 
equations for the composite body has been presented in his 
dissertation. It will be reproduced here for the sake of 
completeness . 


2.3. 1.1 Equation for Points within a Layer 


The first step in finding an equation that will represent 
nodes within a given layer is to substitute equations (20) 
and (21) into equation (11). 


0 i, k, j ~ °i,k, J « 6 l+l,k,j' 2 °i,k.J + °i-l,k,J + °i+l,k,J~ 20 ifk. j 


At 


2L 


(Ax ) 2 


+ °i-l, k, j e t.k+l,j" 2 6 i,k,j + 6 i,k-l,J + 0 i,k+l.j' 2 °i,k,j + 6 i,k-l.J 


(Ay ^) 2 




P J C pJ T ref 


J - 1, .... 5 


(23) 
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, the desired nodal temp- 


This equation is solved for 0. . , 

D 

erature . 


j 

A j { 2 (Lfii) 2 ^° 1+1 > k 'J ” 1-1 , k, j “l+l,k,J T “i-l.k.J 


+ e ,_, w .+ O, , , + °A , J 


+ 2(Lfl ^ } 2 16 i,k+l.J + 8 l.k-l.J + 6 i,k+l.J + Vk-I.J 1 + ®j U i,k,j 


q/'/LAy, 

+ — J-z J— ) 1-2, .... n-1, j = 1 , .... 5 

P J C pJ T rcf 


(24) 


The constants in equation (24) are defined as 

l 


j 


-L+ ■ J + J_ 

At (LAx) 2 (Lfiy^)^ 


(25a) 


J At (Lfii) 2 (Lfiyj) 2 


and 


(25b) 


q” ■ the heat source per unit area, and equals q^LAy^. 


Equation (24) is valid for all the nodes within the compos- 
ite body except at the centerlines of the heater (l = n) and 
the gap (1 = 1). At these points, the insulated boundary 
condition, equation (17), must be satisfied. Finite differ- 
encing equation (17) results m: 


u 


0, k.j 


0 


2,k.J 


J - 1 6 


(26a) 


9 n+l,k,j 8 n-l,k,j 


1, 


(26b) 


Figures 3a and 3b show grid points at the centerline of 
the gap (i = 1) and at the centerline of the heater (i = n), 
respectively. Nodal equations for these points can be de- 
rived by substituting equation (26) into equation (24). 
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At the gap centerline: 


i = x, e 


A. {- 


J 12 0, , + 2 8* 1 + 1 


‘• k 'J J 2(LAx) 2 2,k,J 2,k,:1 l (LAy^ ) 


- >2 ‘“l.k + l.j 


+ o, , , + e 


+ u 


l.k-l.j l,k+l,J l.k-l.j 1 j 1 , k, j 


J - 1 5 


) + B, U, A , + - 1- - - ^ ) 


p C T , 
J PJ ref 


(27a) 


At the heater centerline: 


i = n,0 . = A { J-r [20 + 20 " ] + — 

J 2(I.Ax) 2 l,k,J ’ 2(1 Ay ) 


2 (0 n,k+l,1 


+ 0 +0 A f6,,J + B0.}j = l 5 

n,k-l,j n.k+l.J n.V-l.j J n,k,j 


(27b) 


2. 3. 1.2 Equations at the Layer Interfacial Points 


The interfacial boundary equations are derived by substi- 
tuting equation (22) into equation (15b), and using equation 
(15a). Allowing (i,k) to denote the grid coordinates of an 
interfacial point yields: 


e i,k,j ‘ °i,k,j+i ’ °i,k,j ‘ Yk.j+i * J ‘ 1 5 


(28a) 


( l»k+I,j fik-l,j l.k+l.j f . k-1 ,j ^ a 

4 Ay 

J 

„ _ K ( Yjjtlall 1 * °l,k+l.i+l ~ e i,k-l,j+l ) 

J+1 4 A *J+\ 


(28b) 


It must be noted that temperatures 0^ e i,jc+l,j 
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9 i,k-l,j+l ' anci ^i^k-ljj+l are fictitious and have to be 
eliminated. This is accomplished by writing equation (24) 
for layers j and j+1. 


i,k+l,j + Yk+l.j Aj 0j 1 u i,k,j " A j l 2aA -)i ' U i+»A.J 


2(LAyl 2 


( 0 . 


- A, l- 


(0, 


A A ^ 


"l-l.k.j + Vl.k.j + M't.w,,) + 


--T 


■f 0 


2(LAyj) 2 1,k-1 'J 


A 0?«-Ay 

+ B . 0 , . + -JL J- 


j i.k.J 


P, c .T 
j pj ref 


-]} 


(29) 


+ 6 


2(LAy uiii {e 




-i.k-l.j + l • V i.k-l.j-M “ A J+1 a j+1 lD i,k, J + l - A j+ 1 [ 2(u - )2 


<6 t+l,k,j+l + 0 l-l,k.J+l + °i+l,k,j+l + °l-l ( k. J +l ) 


“j+1 


( 0 . 


2(LAy j+1 ) 2 ' *» k+1 -J +1 °l,k+l,j+l' 


+ B 6 A + 

j+i i.kj+i + - 

J+1 PJ+l ref 


1 ) 


(30) 


Combining equations (28), (29) and (30) yields: 


8 


, » k j . ^ ( . a ^ 2 , 

2 ' , .-.2 .. ,.-.2 1 


i k J 2(LAy i )2 K i+i Ay i 2(LA y 1 +i ) " (Ai> ‘ K j Ay j+i (A *> 

A J “j ^ Ay j + 1 A j + ) « j + l 


l *W.J + Vl.k.J + 0 l+l,k. J + Vl.k.J 1 + 2 l °i,k-l,j + e i\-i,j< 

2 K Ay 2 (LAy ) 2 

+ J- [0 +0 1 + [ J B 

K l.k+l.j+1 l ,k+l ,j + l 1 . J 

K J Ay j+1 J 

2 (UAy . i> 2 K Ay 2 (LAy ) 2 q"/I.Ay 

+ ‘ 1 1 1 B ) H + J 1 

U i+L K . Ay i+1 J+ ’ ’ J "j D C T 

J+1 J j+1 j pj ref 


20-Ay > K Ay q" +1 /LAy 

+ 1 J ; J J J ) 1=2, ... n-1, J=1 , 5 

u j + l P j + 1 C pj+1 T ref 


(31) 
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This equation is applicable for the nodes at the interfaces 
between the layers. Once again, there are fictitious temp- 
eratures that must be eliminated at 1 = 1 and i = n. Sub- 
stituting equation (26) into equation (31) yields.- 
At the gap centerline: 


1-1, 8 


( 2 [- 


(Ay, ) 2 


1,k,J 2 (1-Ay , ) 2 K , Ay 2(LAy (jl .) 2 (Ai) 2 


’ y i+i 1 J )4i' 


J J K j Ay j+i A j+i °j+i 


(Ay 1+1 ) K Ay, 

+ i-= <1 16 ■+ e a 1 + 216 

(Ai) 2 K Ay 2,k,J 2,k,;t 


+ \k-l,j J 


2(Uy,) 


K Ay 

+ 2 -lil — -i- [ 6 + o A ] + [ ^ — b 

- 1 l,k+l.j + l U l,k+l,J + l J 1 a “j 

K j Ay j + 1 J 




>1 K J 1 ’ J Vpj’r.f 

, Vl * ; 1 , 

“j+l K J Ay j+1 C J+1 C pj+1 T ref 


(32a) 


At the heater centerline: 


l , (Ay i } 

i- n , 8 - — j ^ 5 —( 2 ( J— 

2 (t ‘ Ay l ) 2 <L A y l-n) * 1+1 &Y 1 (A *> 

a 1 a A K Ay 

J J J+l A j+1 j y j+l 


, iSttL’ . V> »J, ,0 

(Ax) Kj Ay j+1 


h-1 , k, j + Vl.k.J* + 2 [e n,k-l,j + Vk-l.j 1 


2K 1+ i Ay i a 2<LA V 2(Ld y,+i> K i+i 

+ ] — [0 .. +8, ] + 1 2 B + -^±1 J 

K j Ay j+l " ’ J " >J “ J J “ j + 1 K j Ay j+1 

B , o a + q l' LA ~ y i - + 2(LAy j+i )2 _ K j+i % q "i+i /LAy i +i , } 

J P j C pj T ref “j + l K j Ay j+1 P j + l C pj+l T ref 


(32b) 
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2.3. 1.3 Equations for the Subs t rate- Ambient Interface 

A representation of a node at the inner ambient interface 
is shown m Figure 4a. At the inner surface, j = l. Allowing 
(i,k) to denote an interfacial node at the lower surface, 
the convective boundary condition, equation (16a) becomes: 




°1.0.1 1 °1.2.1 ' 
4 ay 


h bl L ( 


0 +6 

hlxl ii±A . Q > 
2 bl' 


(33) 


Once again fictitious temperatures are encountered, 

0. „ . and 0^ n , , and must be eliminated. Using equati 

1 , 0,1 1 / 0,1 

(24) and letting k = 1, } = 1 and q "= 0, yields: 

D 


on 


+ o . . . . + o . A . . . + p . . ,1 


"l.l.l = A 1 1 2 (LAx) 2_ t ° i+1 * 1 > 1 ' ”1-1. ' '3+1. 1,1 ’ 1-1, 1,1 

* ",.i.i * *i.o.i * *i?i.i * 

2 (LAy^) 


(34) 


Substituting equation (34) into equation (33), 


('iy 1 ) 2 . . 

1 - 2 ( °i+l,l,l + + 0 i+l, 1,1 + °l-l,l,l ) 


1,1,1 2(LAy 1 ) 2 2t. b ,I (Ax) 

" v s r + ~ V “ 

2 (LAy ) 2 2 h LAy 

* 2(0 i,2,i + •£.!> + B 1 * 4: ~ 5 \.1.1 + 


A K bl LAy l . , . . , 

— — °bi } 1 2 ’ n " 1 


(35) 


This equation is correct for all inner surface-ambient in- 
terface nodes except those at the gap and heater center- 
lines. For i = l and i = n, equation (26) must be substi- 


tuted into equation (35). 
At the gap centerline: 
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i-i, e 


1,1,1 Jaiy/ 2h bi LA yi (A *> 

~\ T r+~i 


( ^1 )2 A 

12 hr I 0 2,1.1 + + 


2 [0 1>2(1 + 0 


A , + f 2 (L ^l )2 g 2 h bl L ^l . „ A + 4 \l L % , 

1.2.1* 1 a 1 1 K x 10 l,l.l + ~Kj; V 


(36a) 


At the heater centerline: 


i**n, 0 


t^) 2 

{2 — rV [b_ 


+ e “ , , ] 


11,1,1 2(LAy 1 ) 2 2h bl LAy 1 ' (Ax) 2 11-1,1,1 n -l.l.l’ 

~\ 7 T + h 

. 2(LAy ) 2 2h LAy , 

+ 2 [0 + 0 A . 1 + ( i — B. ^ i] 0 A , 

n,2,l n,2,l a. 1 K. n,l,l 


Ah hi LAy i 

+ _Ai 1 o . 

K ] I>1 > 


(36b) 


2.3. 1.4 Equations for the Abrasion Shi eld -Ambient 
Interface 


When there is no ice on the abrasion shield, convective 
cooling may occur. A representative grid for this condition 
is shown in Figure 4b. Let j = 5 and allow (i,m) to denote 
the interfacial node. Representing the convective cooling 
boundary conditon, equation (16b), in finite difference form 
yields : 


A A 

_ K ( 6 l,nH-l,5 ~ 6 1 , m- 1 , 5 * l.nH-1,5 l^Blbi) 

5 * ay 5 


6 < . + 8, . 

h L ( ■ » m l 5 t»°i 5 _ 0 ) 

"*-•> *• ' •» b2' 


b2 


(37) 
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Temperatures 9i /m -l,5 anci e i ,m-l , 5 ajre fictitious and must 
be eliminated. The procedure is the same as for the inner 
convective surface. Using equations (24), (16b), and (26). 


0 


1 


1,0,5 2(Uy 5 ) 2 2h b2 L4 y 3 (**> 

A 5 °5 " + 


<0y 5 > 2 

{ 2 l °i+l,m,5 + °1-1 ,<n , 5 


+ °l+l,m,5 + °l-l,m,5 ] + 2 l0 l,m-l,5 + 

2(LA; ) 2 2h L«y' 4h b2 LA y5 

+ [ 5_b ^ o* . + i« u ) 1-2. ... n-1 

1 5 i,m,5 b2 


(38) 


At the centerline of the gap: 


i (A y 5 > a 

1 - 1 , o { 2 — V- [0, m , + o ) 

l - m ’ 5 2(LAy 5 ) 2 2\ 2 Lay 5 (Ai> 2,0,5 ’ ’ 

A 5 u 5 K 5 


+ 21U + 0, , .1 + ( 

1,01-1, 5 l,m-l, 5 


2(LAy ) 2 2h LAy 


h2 5, A 


> 5 B '> K 5 ' °1,»,5 


41. L Ay 

+ iT “ 0 b 2 > 


(39a) 


At the centerline of the heater: 


i*n, U 


1 


(4y 5 ) 2 

{ 2 V jo 


+ u 


n,n1,5 2(Ay 5 ) 2 2h fe2 LAy 5 (Ax) 2 n -l.“. 5 n-l.m.5 

A 5 °5 K 5 

2<la 9 ) 2 Zh h2 LA y s 

+ 2(0 , c + 0 ° , C J + [ 2_ B ij o «. 

n,m-l,5 n,n-l,5 a, 5 K, n,m,5 


+ kT ®b2 } 


(40b) 
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2.3.2 Finite Difference Equations for the Ice Layer 
2.3.2. 1 Type 1 Element, Interior Node 


Up to this point, a substitution method has been used ex- 
clusively to derive the governing finite difference equa- 
tions. There is, however, an alternate energy approach that 
uses a nodal energy balance, and has been shown by Arpaci 
[5] to be equivalent to the substitution method. This al- 
ternate method considerably decreases the amount of algebra- 
ic manipulation necessary to reduce an equation to final 
form. It is also easier to physically understand what the 
equation represents. For these reasons, the nodal energy 
balance approach will be used from now on. 


The type 1 element has conductive interfaces on all sides 
as shown in Figure 5a. Writing an energy balance for the 
nodal volume of Figure 7a, yields: 


q 4 m q l + q 2 + q 3 + q STOR 


(41) 


where 


v aT 
ax 


£i 




, 51 

ri 8x 


ri 


ay 6 


a - -K 21 

q 3 ' K ti ay 


ax 


tx 


, - _k 51 

^4 K bi ay 


AX 


bi 


aH 


q STOR = Ft Ax 6y 6 


( 42 ) 
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with the subscript li denoting the left interface, n denot- 
ing the right interface, ti denoting the top interface and 
bi denoting the bottom interface. The thermal conductivi- 
ties at each interface vary with the phase present. Chao 
[ 4 ] approximated this effect by replacing these conductivi- 
ties with average values depending upon the phase of adja- 
cent nodes. The same method is used here. 

= KW = Average conductivity between nodes (i,j) 
and (i-l,j). 

K r ^ = KW 2 = Average conductivity between nodes (i,j) 
and ( 1 + 1 , j ) . 

K ti = KW 3 = Average conductivity between nodes (1,3) 
and ( 1 , j + 1 ) . 


= KW^ = Average conductivity between nodes (i,j) 
and ( i , 3 - 1 ) . 

The partial derivatives m time and space of equation ( 42 ) 
may be approximated using the Crank-Nicolson finite differ- 
ence algorithm at the half time step: 


3T 

ax 


£ 1 


II 

ax 


n 


8T 

ay 


ti 


T i,k,6 

<0 

X. 

1 

H 

f- 

T A 

+ T i 1 k , 6 - 

t a 

J 1 - 1 , k , 6 


2 

AX 


T i+l,k,6 

T i,k,6 

_A 

* J l+l , k , 6 

- t a 

T i.k.6 


2 

AX 


T i ,k+l , 6 

T i ,k,6 

-A 

- T* 
l i,k,6 


2 




+ O(Ax) 


+ 0(ax) 


+ 0Uy 6 ) 


( 44 a) 


( 44 b) 


( 44 c) 
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3 T 

ay 


bi 


^i,k,6 ” 1 i , k- l , 6 + ^l,k,6 ” ^i,k-l,6 
2 ay. 


0( ay 


6) 


( 44d) 


|H . "i.k .6 , out) 2 


( 44e ) 


Substituting equations (43) and (44) into (42) yields: 


Ku ( 1 1 jj ' 6 l- 1 , k , 6 , 

KW l i 2 ax ’ Ay 6 


T A _ r A 

ifW A ( 1 i k i ^ 

* KW i 1 2 U 1 


Ay 6 


‘h 


-kw.,( 


T i+ 1 , k , 6 " T l,k,6 1 ... 

‘ Ay 6 


T*- 

a, *i+l,k,6 i,k,6, ... 

2 m } &y 6 


KW 


KM ” T l,k, 6 v 

" KW 3 ( 2 Ay 6 


AX 


T A , _ 

K ^ 1 ^ 9 i y k , 6 % 

KW 3 C 2 Ay 6 } 


AX 


-KW^( 


T l,k,6 ~ T 1 , k- 1 , 6 
t ay 6 


) AX 


T - T 

KW A (— l - ,k ? 6 2Ay 1,k ~ lf 6 


) ax 


"TO II 


( Hl,k,b at Hl ’ k,6 ) Ax 


Ay 6 


(45) 


Substituting equation (45) into equation (41) and converting 


to dimensionless form yields: 

at T 


e 


h , = h“ . , * r -^~ ( KW,, 

i , k , 6 i,k,6 ” ' 2 


(ay 6 )‘ 


+ KW? (liliillii - 9l,k ’ 6 ) ♦ KW- 

4 (ay 6 ) 2 2 (ax) 2 
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e a e A 

+ KW A ( 1 - 1 > k > 6 i . ’ . k ’ 6 ) * kw 

<^6> 

e a - e A 

♦ KW* ( - j -JSli-l 6 _ J _L. i l.: . 6 + K w 

3 ^y 6 > 


e A _ e a 

+ KW A ( i*~ 1 1 k , 6 1 1 k t 6 

1 (AX ) 2 

KW, + KW. 

+ _1] e, ,) 

( a y 6 ) 2 ,k ’ 6 


l,k+l,6 
3 <**6 )2 


0 i- 1 , k , 6 

1 (A-X) 2 


) - [ 


KW X ♦ KW 2 
(a *) 2 


(46) 


An equation for each phase is developed by incorporating 

equation (14) into equation (46) and solving for H. , 

1 , K , b 

For the solid phase : 


Vk ,6 " (1 


At 


* P u — P 

2 p s c ps L Ux) 


KW, ♦ KW, KW, + KW„ , 

c — - — - — 2 r 1 


(Ay 6 ) 2 


At T 


,H A , , + [KWjj «. Kwjd^.6 “ ^ k ' 6 

’ ,t> 2 L 2 ( Ayj 2 " ( Ay 6 ) 2 


+ KW 1 * 1 , - k i 6 + KW* (— — 1 : k i 6 - - 1 ’ k 6 )+ KW, - J ■ k * . -L ’ 6 

( a 5) 2 2 ( A *) 2 3 ( Ay 6 ) 2 

+ kw* ( -l-i — * 1 ■ 6 — - 1 : k ’- 6 + KW, 9l - 1,k ^ ♦ kw* ( 6 l ~ 1 -t k ^ e ,h; k ’ 6 ) ] ) 


(Ay 6 > 


( AX ) 


(AX)‘ 


(47a) 


For the melting point : 

At T 


H 


i,k,6 - i , k , 6 T 


(KW, 


<*'V r 


+ Kwt ( 9 i -t k . - 1 ^ 1 ' k -' 6 ) ♦ KW, - — ]_ 8 - y 1 - 6 


(Ay 6 ) 


2 


♦ KW* ( Il f liiL ti ♦ KW, ■ * f k - | ~ 1 i 6 

2 (Ay 6 ) 2 3 (Ay 6 ) 2 


a e 


KW. 


f °l,k+l,6 ~ e l,k,6 ) + KW ”1-1. k, 6 

/ . * v 2 1 / . " » 2 


(Ay 6 ) 


( AX } 


33 



(47b) 


A e f i u ft - e i u R KW. ♦ KW- KW- + KW, 

+ KW ^ ( l,k,6 _ t_l 2 + _J 1 3 e ) 

1 (ax ) 2 (ax ) 2 (Ay 6 )^ 


For the liquid phase : 


H i . k , 6 ' 


At 


2 P £ C p* L 


KW X + KW 2 KW 3 ♦ KWjj _ x 


(ax)‘ 


(Ay 6 ) 2 


{H i k 6 * [KW y -- k " 1 ; 6 ♦ Kwf( 1 k- 1 f 6 jLlJLiJL) 

l.k.O 2 L d 4 < a ? 6 ) 2 4 ' 2 


(flye)* 


KW 1>1 » k . ’ 6 + KW^( - - ■ 1 -j k , 6 f A P .k,6 ) 

2 (a5) 2 * (a*) 2 


+ KW. - 1 - ’ k * 1 \ -- + KW*( I i k->,1 I 6 - JL’JL’1 ) a- KW. *~Aj k >. ( 

3 (ax ) 2 3 ( Ay g ) 2 1 (ax ) 2 

+ KW,( - - '.-» k ’ - 6 — ^ Ls JL>A ) 

1 ( AX ) 2 


KW. •* KW- KW- ♦ KW. . . L _ 

[— p- 5 * - 3 ". 2 3 (C Ds 0 m ♦ J_ + el]) 

■ = ' 2 (aP 6 ) 2 c p£ PS m VeT m 


(ax)‘ 


(47c) 


These equations are identical to those developed by Chao for 
a purely conductive node, and are second order accurate. 
This method is equivalent to the substitution method, but is 
more direct and easier to apply. 

The insulated boundary conditions at the heater and gap 
centerlines are handled within the computer algorithm. For 
example, at the gap centerline (i = 1) the boundary condi- 
tion is expressed m finite difference form as: 


0 »,k,6 _ °0,k,6 
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The fictitious quantity 6 n , , is replaced by 8. , ,, in 

U f K ^ u c. / K f O 

equation (47). At the heater centerline (1 = n), the bound- 
ary condition becomes: 

0 n-l,k,6 “ °n+l,k,6 

The nonexistent quantity, 0n+l,k,6, in equation (47) is re- 
placed by 0n-l,k,6 • This substitution process is handled 
entirely within the computer program, and reduces the total 
number of equations. It also simplifies the numerical model 
considerably . 


2. 3. 2. 2 Type 2 Element, Ice-Ambient Interface 

This type of element is shown in Figure 5b. The same en- 
ergy balance still applies, but one surface is now convec- 
tive. The terms of equation (41) become: 


q i 


q 2 


q 3 


KW T l-l,kj6) 

KW 1 ' 2 ax 

K u a ( 1 i k ^ 1-1 , k , 6 ^ 

KW 1 1 2 lx 
_ KW / T l+l. k )6 

rw 2 v 2 &x 
A _ 

i/ijA , 1+ 1 1 k , 6 i , k , 6 v 

KW 2' 2 AX 

, T i , k , 6 + T 1 i k , 6 

h b2 ( 5 


* y 6 


Ay 6 


Ay 6 


ay 6 


X b2 ) 


AX 
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( 48 ) 




-KW U ( Tl ' k '% AX 

d a y 6 


KW ^ 1 ^ 6 Tf t k ~M ) AX 
4 2Ay fi 

, H i,k,6 " H i,k,6, 

q STOR “ ( !— SI L ~‘~ * Ax iy 6 

Substituting these into equation (41) and combining the re- 
sult with equation (14) produces the final result for each 
phase . 

For the solid phase : 


Vk,6 ■ 


At r KW 4 KW l * Kw 2 h b2 L ,,.i 

2 1 L r~77 + — 2 + 3) 


2L "s c ps Uy 6 } 


(AX)‘ 


Ay ( 


(H 


A . * re ^ neu l , k- 1 , 6 tu A , e l , k-1 , 6 **i,k,6* 

‘•“■ 6 2 i. 2 ‘‘ 7 ^^ ^ (.; 6 ) 2 


e, 


KW ■ ■- 4* x i . k ^ - KWp( 8 ~ ^ f k t" 6 ° 1,k ' 6 ) ♦ KW, e i-l,k,6 


- 0; 


a e. A , . , - ef , , h KO l 

♦ KW. ( 1 ~ 1 4 k . »^ ... fjL . k /.§ ) + _Jlf! 

( AX ) 2 


(2 0b 2 - e ? k 6 )]) 

A y t 1 f K , O 


(49a) 


For the melting point : 

H . „a * At T ref ,„ w e l,k+l,6 . „ W A , 9 l,k-l,6 ~ e i,k,6 , 

i,k,0 i,k,b 2l 2 Uy 6 ) 2 q Uy 6 ) 2 


KW- k l - ♦ KW, ( 9 l->, 1 > k i - ^ -i k ’ 6 ) 

2 (A *) 2 2 (AX ) 2 


KW ei ' 1 ,k j , 6 + KW* (lizLJili L L k , 6 ) 

1 (Ax ) 2 1 (AX ) 2 


e? , .. c - 0? 


h.-L . KWi, 

— — — — ( 2 e b 2 - 0* fi)- [ ? 

Ay 6 b2 i>k ’ 6 (Ay 6 ) 2 


KW X + KW 2 h fc2 L 


(a5) 2 


* y 6 


1 

m 


(49b) 


For the liquid phase-. 
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"i.k. 6 ■ 11 ♦ rr 


At 


KW„ 


2L- Pf C pt (ay 6 > 2 


KW i + KW 2 h b2 L ,.-l 

♦ 31 , 


(ax) 


Ay# 


At T 


££l [KWi, 0i «. KW f; ( ° 1 i k ~ 1 i 6 — Z 1 ! - ,6 ) 


(H 


i,k,6 + 


2L 


* ( ^6 )2 


(Ay 6 ) 2 


♦ KW + KWp ( e i-l , k , 6 "M’^ ) + KW ei - 1 ’ k ^ § 

2 (AX) 2 2 (AX) 2 1 (a5) 2 


V1 , a , °i- 1 , k , 6 “ e i,k,6 x . "b2“ „ Q 

+ KW i ( L ~ -x ~2 } (2 b2 " 9 i ,k , 6) 


hf -.1- 


( AX ) 




KW 4 KW, + KW 2 h b2 L 

[ -7-T 2 + - — r - + — : 1 [ T“ (C PS °m * r-f 


+ JlL.) . e _]]} (49c) 


1 (>' 4y 6 C pt 


(Ay 6 ) (ax) 

The insulated boundary conditions at the heater and gap cen- 
terlines are handled within the computer algorithm in the 
same manner as the type 1 element. 


2. 3. 2. 3 Type 3 Element, Ice-Ambient Interface 


This element has a convective surface on the left nodal 
interface and conductive surfaces on the others, as shown 
in Fig. 5c. The energy terms of equation (41) become: 


q i 

q 2 

q 3 


♦ T, 


h b2 ( - > k ' 6 „ 1,k ’ 6 - T ko ) 


b2 ■ 


&y 6 


-KW-( Iul > k >| A , Tl,k -^) 


2 Ax 

K w ,( 2 r- ax ' &y 6 


Ay 6 


-Kw / 1,k * 1 '; Tl,1<L -) AX 

3 Z Ay fi 


KW 


— “r 

a, ' 1 , k+ 1 , 6 1 , k , 6 ^ 

2 Ay, 


AX 
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(50) 


-KW { * > k ’ b » ^-1 , 6 ^ 

K *V 2 Ay 6 ' 


KW* ( — j-t-iL t_ §. , ) Ax 


( Mlfkf6 At Hi - ^’ 6 ) AX Ay fi 


q ST0R 


Carrying out the substitution of these terms into equation 

(41), incorporating equation (14), and solving for 

H . „ yields the equations for each phase, 

i/k, 6 

For the solid phase : 


Vk.6 - 11 * 


At 


KW- t- KWj. KW- hb-L . 

C—2 ^ ^ ♦ — S-]}" 1 • 


"5 L - 5 

2 p C L (Ay,) ( x) c a x 

s ps o 


{H* , + [KWj, 1,k ~ 1 ’ 6 + Kw£ (-iiitllti ’ 6 ) 

1>k ' 6 2L 2 4 (Ay b ) 2 Uy 6 ) 2 


* KW, 0:l > k t^ 6 ♦ KW* ( 0 ±!. k li.JL b f 1 i f k i 6 ) _ KW- lilitii- 6 

3 / A - ,2 2 u5) <? 


3 <*y 6 )2 


(Ay,) 


♦ KW* ( 0 1^.X.6 ~ 9 ".^) + ^2 L 
2 (Ax) 2 AX 


T" ( 2 °b2 " 0 ?,k,6 )]) 


(51a) 


For the melting point ; 


H 


i,k,6 


At T 


H 


i > k , 6 


ref 


6 


KW 


2L 


1 , k-1 , 6 
6> 


* <*y fi ) 2 


+ KW* (.. A; . * ! 1 ) 6 i|k - ’^ ) + KW- 9 ^ -- 4 - 

(*y 6 ) 2 3 (Ay 6 ) 2 

+ KW* (■ ki k *h . ^ l . ’ . k - ’ 6 + KW- ■ 1 — i-j 5 - ’ . 6 

3 ( y 6 ) 2 2 ( X) 2 


KWt 


- [ 


KW 3 + KW 4 


+ ^ (2 e b2 - e* k 6 ) 

(AX) 2 *5 b2 1 ’ k ’ 6 

KW 


hb 2 L 

AX 


- 5 + - o * 


(Ay,) 


(a*)‘ 


hb 2 L 

AX 


] 6 m f 
tn 


(51b) 
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For the liquid phase : 


H i , k , 6 = 11 * 


At 


? L f 

2p. l Uy 6 r 


r KW 3 * KW| i KW 2 h b2 L 11 -l 
L 5 - + =5 ♦ — J ) • 


l pi 

At t e 


(a5) 2 


AX 




(H“ + i— [KWj. - - t— ♦ kw£ ( 6 L ti i i i ) 

’ 2L 2 4 (Ay 6 ) 2 14 ( Ay 6 ) 2 


+ KW, - + KW* (JjJStlti 1 » k » 6 ) + KW, ■ l U 1 i %> 6 

J (Ay 6 ) 2 3 (Ay 6 ) 2 2 (AX ) 2 

* KWp (hlLtJiA J±l±l±) + -4- (2 9b2 - e i,k,6 ) 

2 (AX) 2 ax 


KW 3 + KW^ 


(Ay & ) 2 (ax) 


KW 2 h b2 L , r 1 * L r 

5 + j L- (C 0_ ♦ TK 

•■' 2 AX C , PS m T r ef 
pi 


) - 0 ] ) ) 
m 


(51c) 


The insulated boundary condition at the heater centerline is 
handled within the computer algorithm in the same manner as 
the type 1 element. This type of element can never occur 
over the centerline of the gap ( i=l ) . 


2. 3. 2. 4 Type 4 Element, Ice-Ambient Interface 


Figure 5d shows a type 4 element. It has a convective 
surface on the right nodal interface and conductive surfaces 
on the others. The terms of the energy equation are: 

- T, 


» K« X ( i. k .6 ~ x Ay 6 

«■ K W* ( Tl,k 4 aI 1 ' 1 - 1 — > A *6 


hb 2 ( 


T l,k,6 * T ^. k L,6 _ ) Ay 6 
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* y 6 


Kwt( Tl ’ k ' flf | AX 


T Ay 6 

g - ( 1 > k ! § 1 , k-l , 6 ^ 

q 4 " K *V 2 Ay, ' 


AX 


I, yA / : . k , 6 i , k - 1 , 6 

W 4 l PAy 


6 - H 1 , k , 6 x 

q STOR 8 ( It ’ * x * y 


A 

±j 

6 

A 


) AX 


(52) 


Replacing the energy terms in equation (41) with equation 
(52), incorporating equation (14), and solving for H. , ,- , 

1 , K , O 

results in the governing field equations for each phase. 
For the solid phase : 

At 


H i , k , 6 8 * 


•5 [■ f- + m ■* 

2 p s C ps L <* y 6> (AX> 4 5 


KW. ♦ KW, KW, hb.L , 

^ 1 - 1 - 2 ]» , 


(H 


At T 


l,k,6 + — 


ref 


e 


[KW„ 


t^ k-l . 6 + kh a (. ^l t k-l F 6 v l t k l 6 ^ 

( avc ) 2 


( Ay 


9 ' ’■ 1.6 ~ e : 
( a? 6 )J 


KW ♦ KW* ( l -t. ** 1 ! 6 " ■ 1 r k , r 6 ) + KW. WtM 

3 (Ay fi ) 2 3 (Ay,) 2 1 (ax) 2 


♦ KW* ( ” 1 - 1 » k i 6 lixiLtl) + , ht>2L 


(2 - 9? „ ,)]) 


( Ax ) c ax »>2 1 ,k,6 


For the melting point : 

At T 


ii ..A . ref fwu i.k-1.6 

H i,k,6 - H i i k , 6 + 91 S~- {KW 4 


* KW* — *i k i 6 ) + KW 3 


2L 

A 


( ay 6 )2 


(Aye) 


KW* (- f ’ k * 1 l 6 , 1 F k l- 6 .) ♦ KW. 1 

3 (A?,) 2 1 (AX) 2 


(53a) 
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* KW 1 ^ (e b2 - e^ (kt6 ) 


(a*)‘ 


- C- 


KW, ♦ KW. 

< a_ v 2 


hb 2 L 


KW X 

— r~? + 

( AX ) 2 AX 


AX 


^ 0 m> 

u) 


Fo r the liqui d phase : 


(53b) 


H i , k , 6 » 11 * 


At 


KW, + KW 3 KW, 


2 p i C p* ( ^6 )2 (A *> 


"•"1 " bS ^,-! 

♦ — —* + — — D > , 

(ax) ax 


— pL CKW, 1 > k - 1 t 6 ♦ Kwf ( 1 > k - 1 > 6 : 
2L 2 (Ayg) 2 Uyg) 2 


KW 


_ 9 l J k, l 6 + kw A ( 9 ?,k+l,6 - ^k f 6 } + ? 1 -. 14,6 

J ' 2 * ( Ayg ) 2 ( ax ) 2 


< A V‘ 


+ kw a ( -l-l f k f 6 • < i J k t t ) , ^ (2 . - e, k 6 ) 

1 ( a *) 2 ax ^ 2 1 > k > 6 

r KW, + KW 3 KWj _ L f 

(Ay 6 ) 2 (AX ) 2 Ax C pS ref 




) - e m ]]> 

m 


(53c) 


The insulated boundary condition over the heater gap center- 
line is handled the same as with a type 1 element. This 
type of element can never occur over the heater centerline 
(1 = n). 

2. 3. 2. 5 Type 5 Element, Ice-Ambient Interface 


This element has convective surfaces on the top and left 
side interfaces as shown in Figure 5e. The energy terms be- 
come : 


q l = h b2 


( Tl,k ' 6 l — t . 1 r k i 6 - T k . ) 


b 2> 4y 6 
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-KW 2 ( 


T i+1 , k . 6 ' T i,k,6 


) *y 6 


KW^f i^l » k , 6 i,k f 6v 

KW C ( TTi } Ay 6 


. , i , k , 6 i , k , 6 « x 

d2 ( 5 T b 2 } Ax 


Qij = -KW^ ( T± 8 - { - ' 6 j A y 1,k ~ 1 ' 6 ) AX 


A A 

KW T i , k- 1 f 6 


Q _ ( 1 * k , 6 i * k , 6 x 

q STOR ~ ( At ' Ax Ay 6 


(54) 


Placing these terms into the energy balance, equation (41), 
incorporating equation (14), and solving for ^ 6 yields 
the governing field equation for each phase. 

For the solid phase: 


Vk,6 = {1 * — 


At r ^2 ,11 ni -l 

, — [ 5 + S + h. -L <— + » 

2 ft C ps < Ay 6 } AX Ay 6 


fH A , + ^ [KW, 1 1 k ~U - 6 + KW? ( 1 f k ~ L- ^ ~ b ) 

1 - k - 6 2L 2 4 (Ay 6 ) 2 4 (Ay 6 ) 2 

e e A - e A 

KW_ + kw a — _Ll 

2 (ax) 2 2 (a*) 2 


+ h. L(2 e. , - e A ,) (— + — )]) 
52 b2 1 - k - 6 ax Ay 6 


(55a) 


For the melting point: 


** 1 , k , 6 * ^ 1 , k , 6 + 


i.k-1,6 


+ Kw£ (1Ll!<zLiJ — i,k,6 ) + 0 i+l,k / 6 

(Ay*) 2 2 (ax) 2 
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* KWp ( i - 1 ’ k - ’ 6 6 ) ♦ h L (2 


’? ' 5 — 1 — * n. , l - e. . ,) 

d (ax) 2 b2 b2 i.K.o 


1 x r KW <, KW 2 

- [ ♦ ^ p 

AX Ay 6 (Ay 6 ) c (Ax)* 


(— 


h b2 L (J r + -=-” V 

ax Ay 6 


(55b) 


Fo r the liqui d p hase : 


H i , k , 6 = (1 * 


At f KW « KW 2 


"5 ♦ ~3 * h h2 L ^~ + - 1 ■] )~ ] * 

2 p ( C pt L c (aP 6 ) ( a5) 2 ax y 6 


(H 


+ [KWj, 6 1 ’ k ~ 1 - 1 6 + KWy (liiiitlii g 1 - ^ * 6 ) 

* (M 6 ) 2 4 (Ay 6 ) 2 


1 - k ’ 6 ’ 2l 2 


♦ KWp «. KWp < 9 ^>. k » 6 '/ I '*’ 6 ) 

2 (ax) 2 2 (A*) 2 


1 ^ 1 


KW, 


* h b2 1(2 e b2 ■ < k * C -• -» 

b2 b2 1,k ’ b Ax Ay, (Ay,) 2 


KW, 


* . -% * h b2 L( ~ T^ _)]C T~ ( S>s 9 m + T77 ) ' e B ]:,} 
(Ax) ax Ay, C . y ref 

b pt 


(55c) 


The insulated boundary condition at the heater centerline is 
handled m the same manner as for a type 1 element. This 
element cannot exist over the gap centerline (l = 1). 

2. 3. 2. 6 Type 6 Element, Ice-Ambient Interface 


Figure 5f shows the type 6 element which has convective 
surfaces on the top and right side nodal interfaces. The 
energy terms for this situation are: 


T 1 , k , 6 ~ T i- 1 , k , 6 


q l = KW 1 ( — '"‘2 a'x 7 af 6 


) Ay, 


+ KW? ( Tl > - ’^ — -litiJiil) 

1 d aX 


A *6 
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hb, ( 1,k ’ 6 —ZULtl 


T b 2 * Ay 6 


. , T l,k,6 + T f , k , 6 T . 

q 3 = h b 2 ( 5 T b2 ) Ax 


T i,k,6 ' X x.k-1,6, ... 

q H - - KV V TTT^ } AX 

A A 

t 1 ^ » 6 l.k-1,6 . 


, H i , k , 6 - H i.k,6, 

( > Ax Ay 6 


(56) 


Substituting equation (56) into equation (41), incorporating 

equation (14), and solving for H . .yields the three gov- 

1 , ]c , o 

erning equations corresponding to each phase. 

For the solid phase: 


H x,k,6 « (1 + 


2 Ps C ps L Uy 6 ) (AX) 


KWj. KW. i,i 

♦ Kf + h. L(— ■♦■—)] ) _1 , 

* 7. / .3)2 02 - . - 


Ax Ay 6 


At T 0 , . . , . e? , , c - 0 ? , c 

[KW, 1 ' k ' 1 ^ 6 ♦ Kw£ ( -LlJSzLll ■ 1 » k -t- 6 ] 

2L 2 4 (Ay g ) 2 (Ay) 2 


)“ + pi- [KW, ♦ KW“ < . ■ x ’^ u ) 

1,K,0 Pf « 4 (AV ^ 4 (Ay)^ 


KW + Kwf ( 9 l -~ . 1 i t b§ . 9 l i k . » 6 ) 

1 (ax ) 2 1 (ax ) 2 


+ h h _ L ( 2 0 - 0 A ,) (— — )]) 

b2 b2 i>k ’ 6 a5 a? 6 


(57a) 


For the melting-point : 


H i,k,6 " H i , k , 6 + 


, k- 1 , 6 


(Ay 6 )‘ 


♦ kh A ( e i,k-l,6 - e l, *j±) + KW. 

H (Ay ,) 2 1 (ax ) 2 
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e* , . , - ef 


+ KW^ (-— ~ 1 ' k i 6 1 1 k ’ j ) + 

(AX) 2 


h b2 L<2e b2 “ e l,k,6 ) 


i± 


-)- [- 




KW 


Ax Ay & (Ay 6 ) 2 (Ax) 


1 1 * h b 2 L ( ~r * “^ ):l e m ) 
ax Ay & 


(57b) 


For the liquid phase : 


H i >k( 6 - * 


At 


KW, 


KW 


2 p< L " 


* h b2 L + — )]) _1 
d ax Ay^ 


(H, 


KW, 


b2 


At 

,6 + — 

^ref 

2L 2 

[kw 4 

e i,k 
( Ay 


6 i-l,k 

tl + 

* e 

KW fl 

A 

1-1, 

k , 6 

( AX ) 2 




(a5) 

L(20 b 2 

- 6 i 

,k,6^ ( 

1 

- — 

-J-) 




AX 

* y 6 

L(JL ♦ 

— )3 

[JL 

(C 

e - 

AX 

Ay 

c p« 

ps 

m 


, l,k-l t 6 _ 1 

(Ay 6 ) 2 


KW, 


KW, 


ref 


) - e ] ] ) 

m 


(57c) 


The insulated boundary condition over the heater gap center- 
line is handled the same as for a type 1 element. This ele- 
ment cannot occur over the heater centerline (l = n) . 


2. 3. 2. 7 Type 7 Element, Ice-Ambient Interface 


This type of element has convective surfaces on the left 
and right nodal interfaces and conductive surfaces on the 
top and bottom interfaces as shown in Figure 6a. The energy 
terms are: 


= h b2 ( i .>. k -u 6 . r 


1 ^ _ T 


b2 


) a y 6 
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»2 


q 3 = -KW 3 ( Tl > k * 1 >| Ay ^ 1,k ’ 6 ) AX 

- KW 3 (- T - ,k '" - , - g - Ay ^ 1>k ' 6 ) AX 


Q _ -KW ( 1 ' — 1 ^ i , k- 1 1 6 v 

q 4 “ KW 4 ( 2 Ay 6 ' Ax 


K^(--- |k i 6 . T i ,k-l , 6 


O - f Vk,6 ~ H l.k,6 

q STOR * ( aI ) Ax Ay 6 


Substituting equation (58) into the energy balance, equation 

(41), incorporating equation (14), and solving for 

H . , produces the governing equations for the three phas- 
l , k , o 

es . 

For the solid phase : 


H l,k,6 = {1 + 


2 P 3 C P3 L U *6> 


KW, + KW„ 2 h K , L , 

r 3 h b 2 ,,-1 

L ; = — + ; J) 

(Ay ct Ax 


a At T „ 2 h. _ L 9. . . , 

fuA ref r b2 , i,k+l,6 

,H i k 6 5 — L ' 2 e h? - 9 < i, f,' * KW a J 2 — 

i.K.b 2L 2 t - b2 i,k,b 3 (A | } 2 

+ kw a ( ’ 6 ^i- ’ . k ’§ ) + KWj. ■-^■> k ~ 1 ; 6 

3 (Ay fi ) 2 (Ay fi ) 2 


+ KW^ ( 1 ' k ~ 1 f 6 LJLlI)]} 


(Ay 6 )‘ 


For the melting point: 
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h . T ref , 2 h b2 L , „ „a 

H i,k,6 - H i ,k, 6 + ~"2 { (2 e b2 ' e i,k,6 ) 


KW . kw« ( f j jLtUl ' °1 ,K,6 , , KU 


<•»,.* - u»,,‘ 


A 9 1 k-1 6 _ °1 k f, 2 h K5 KW 5 ♦ KW 2l 

- KW^ ( - Ll X-1,6 1,*'* ) _ c b2 + _J Jij Q 

Uy 6 ) ax (ay 6 > 2 m 


(59b) 


For the liquid phase : 


at r 2 h b2 KW 3 * K V,-1 

H = ll «■ P [ ~ + I 5-J* 

2 ft C p* L AX ( * Y * ) 


at T r 2 hb-L . e, . . , 

(H ! k 6 + 5^ C -T < 2 e b2 - 9 i k 6> + KW 3 • ^ 

1,K, ° 2L^ tx b2 i,K,t> i (ay 6 )2 


♦ kw? (_ L iUlt i — e ^- k - » - 6 ) + KW }Y~ X J > 
5 Uy,) 2 


A 0 t k 1 fi - k fi 2 h b? L KH T + KW U 

♦ KW A ( i t k-1 ^6 i i k i s ) + t b2_ + — 3 *•] 

(ay 6 )^ ax (ay^)^ 


^ { ' c p3 0 m * T~^~ ) - 0 n ^> 

C P* ref 


(59c) 


The insulated boundary condition does not apply here because 
this type of element can never occur over the heater or gap 
centerlines (i = 1 or n) . 


2. 3. 2. 8 Type 8 Element, Ice-Ambient Interface 

Figure 6b shows the type 8 element, which has convective 
surfaces on the right, left, and top sides and a conductive 
nodal interface on the bottom. The energy terms for these 
conditions are: 
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, T i.k.6 + T i.k,6 T ^ Au . 

<>! * h b2 ( '' 2 T b2 ) y 6 


q, = hb p ( 


T i.k,6 * T l,k.,6 


- T b2 > Ay 6 


„ u , 1 , k , 6 1 f k , 6 t? \ , u 

q 3 = h b2 ( 1 i -~ g - _ T b2 ) Ax 


a - -KW ( * ’ il ’ ^ 1 1 k- 1 , 6 . 

q 4 - kw 4' 2 Ay & ' 

A A 

vu A r - 1 < ^ , 6 1 , k-1 , 6 , 


q STOR " < Hl,k,6 A ' t Hl ' jk,6 ) Ax Ay 6 


(60) 


Introducing these terms into the energy balance, equation 
(41), incorporating the temperature enthalpy relationship, 
equation (14), and solving for fh ^ g yields the governing 
equation for each phase. 

For the solid phase : 


Vk,6 ■ 11 


At r ** u , , 2 . 

— * 5 L ^ ~L v— ♦ )J) « 

2 Ps c ps L < A V Ax Ay 6 


(H“ c * 


i,k,6 


1 , k- 1 , 6 ^ „ u a , l.K-1 ,6 ~ 8 1 , k, 6 , 


Uy b Y 


+ KW^ ( 


< A h y 


h h ,L (2 e h - - 8. .,)(-£. ♦ — )3J 
b2 b2 ilk ' 6 ax Ay 6 


For the melting point: 
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At T 


H i , k , 6 = H i , k , 6 + 


ref 


KW 


*1 , k-1 , 6 

“ (ly 6 ) 2 


+ KW* j ^ : vr — ? + h b2 L < 2 e b2 - ®1 ,k,6 ) + 

(Ay 6 )‘ : ’ ’ A X Ay ft 


- [- 


KW„ 




(Ay 6 ) 


'fc 2 L ijz * V 


Ax Ay 


(61b) 


For the liquid phase : 


Vk ,6 ■ {1 


At 


KW, 


2 + h b2 


2 Pz L <*V 


h. _ L (~ ♦ , 


ax Ay, 


At T 


£ [KWjj 6l, H' 1 ^ 6 ♦ KW* ( 8l > k ~ 1 > 6 _ 8l -t ^ 


/ ret 

i,k, 6 + 2L 2 - (Ay fi ) c - (Ay 6 )‘ 


KW 

+ ^b? ^ ^ ^ 8. p - Bj . ,) (— ♦ — — ) + [- — * 

b2 b2 ilk * 6 Ax Ay 6 (Ay 6 ) 2 


♦ h.pL(— ♦ — )] [J- ( c 8 - - 1 -) - 8 ] 1 } 

62 « *?e Si p ^ 


(61c) 


The insulated boundary conditions at the heater and gap cen- 
terlines are of no concern because this type of element can 
never occur at these positions. 


2. 3. 2. 9 Type 9 Element, Abrasion Shield- Ice Interface 

The final type of element occurs at the abrasion shield- 
ice interface, as shown in Figure 6c. Figure 7b shows this 
type of node and the associated energies passing into and 
out of the element. The energy balance equation is slightly 
more involved: 
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a 4 « q a + q* + q a + q* + <1 + <1 

q ^ m 3 t ST0R 


The energy terms of this equation are: 


.a „ , T i,k,5 ' T i-1 , k , 5 + T i,k,5 " T i-l,k,5, Ay 5 
q l ’ K 5 ( ’ ~T~ 


i / T i i k , 6 ~ , k , 6 , 4y 6 

1 - KW 1 ( Sax ) ~T 




(62) 


a / T i+1 , k , 5 " T i.k,5 * T i*l.k.5 " T i.k,5 ) Ay 5 

q 2 - ~ K 5 1 2Ax ~T~ 


q* . _Ky ( i-»-l > k,6 l,k,6 > y j 

q 2 - ™2 y 2 ax ‘ 


- KW“ (• 


r i+l,k,6 " T i,k,6, Ay 6 


q z -KW. Tl i k > 6 ) ax - Kwt — Tl ' k ’ 6 ) Ax 


<ln = "Kc ( 


T 1 , k , 5 ~ T 1 , k-1 , 5 *’ T l,k,5 ~ T x ,k-l , 5, 


, H i , k , 6 ‘ H i,k,6, Ax Ay 6 „ , T i,k,6 ~ T i,k,6, AxAy b 

q STOR = ( M At L ^“ ) — 2 + p 5 C p5 ‘ ‘Tt L “ J - ) — 2— 


(63) 


Substituting these terms into equation (62), incorporating 

the enthalpy- temperature relationship, equation (14), and 

solving for H . yields the final governing equations for 
1,JC, 6 

the three phases. 

For the solid phase: 


„ r , P 5 C P5 At , 2 k 5 2 kw 3 

n . c - l J- ♦ * + * 5 L ♦ 5 

p s c ps 2p s C ps L Ay 5 Ay 6 (Ay 6> 
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1 " c 

<aS) 2 
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(64b) 


For the liquid phase : 
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These equations are identical to those developed by Chao (4) 
for the abrasion shield-ice interface, further validating 
the energy balance approach. 
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2.4 NUMERICAL SOLUTION BY COMPUTER IMPLEMENTATION 
2.4.1 Gauss-Seidel Point Iterative Method 

The Crank-Nicolson finite differencing formulation re- 
sults in a set of simultaneous algebraic equations that must 
be solved at each point in the spatial grid for each time 
step. There are several methods available for solving the 
matrices resulting from the set of simultaneous equations. 
Some use matrix inversions, and others use iterative pro- 
cesses. An iterative process is necessary for the ice layer 
because the nodal phase condition, solid, liquid, or ice, is 
not known m advance. The Gauss-Seidel point iterative 
method was chosen because it was already incorporated by 
Chao [4] in solving for the composite body nodal tempera- 
tures. He had also used it for the ice layer. 

The Gauss-Seidel iterative process starts for each time 
step with an initial approximation derived from the previous 
time step. The equations are solved for the desired nodal 
quantity, enthalpy or temperature, as can be seen from the 
final equations for the ice layer and the composite body. 
The initial approximation is improved by passing through the 
spatial grid and solving the equations at each point. This 
process is repeated many times, eventually resulting in a 
converged set of values satisfying the set of equations for 
the spatial grid of nodes. 
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For this study, the iteration process was considered to 
be converged when the difference between the previous and 
present iteration was less than 0.01%. Once this criterion 
was met, the iterative process was started for the next time 
step, using the final converged matrix of values from the 
previous time step as the initial values. 

The Crank-Nicolson formulation results m an uncondition- 
ally stable set of equations. Reference [9] contains a more 
complete discussion of the suitability and stability of the 
method . 

2.4.2 Numerical Program Algorithm 

A complete flow chart for the computer program is shown 
in Figures 8 and 9. It is identical to Chao's since this 
study was based on his work. The mam program is almost 
identical to Chao's, except for a few changes necessary to 
allow incorporation of the variable ice thickness subrou- 
tine. Parts of the original program pertaining to the ice 
layer were replaced by the new variable ice thickness algor- 
ithm. 

All of the features of the original Chao program are 
still present except for ice shedding. Information pertain- 
ing to this has been eliminated from the input data set. 
However, the original shedding algorithm has not been delet- 
ed from the program and could be made operational. 
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Chapter 3 

DISCUSSION OF RESULTS 


3.1 VERIFICATION OF THE VARIABLE ICE THICKNESS ALGORITHM 

The first step m validating the variable ice thickness 
algorithm was to compare it to Marano's [3| and Chao's [4] 
numerical models for the case of two-dimensional conditions 
with uniform ice thickness. These comparisons were made in 
order to make sure that the program correctly handled the 
base case of constant thickness ice. Heater gaps were in- 
cluded in many of the two-dimensional runs. For all cases, 
the variable ice thickness model simplified correctly and 
yielded virtually identical results. 

Another approach was used in validating the results for 
variable thickness ice layers. Symmetrical ice shapes were 
run to see if symmetrical temperature distributions result- 
ed. In all of the cases studied, a symmetrical temperature 
distribtion was observed. 
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3.2 VARIABLE ICE THICKNESS PROGRAM RESULTS 


As a preliminary step towards implementing the two-dimen- 
sional code for a variable thickness ice layer, Marano's 
one-dimensional code was run for a deicer model with various 
ice layer thicknesses. The object of these runs was to aid 
in the determination of relevant changes in thermal tran- 
sients due to ice layer thickness. This information can be 
provided using the one-dimensional model with much less cpu 
time than with a two-dimensional code. 

Figure 10 is the plot obtained from results generated 

with Marano's code using the deicer model of Table 1. The 

temperature rise of the ice-abrasion shield interface is 

2 

plotted against time for a heater power density of 25 W/m . 

and for various uniform ice layer thicknesses. The external 

2 

convection coefficient is 200 BTU/hr-ft -° F. Temperature 
rise is defined as the total rise to 32 °F from the initial 
temperature in the ice layer. Phase change is not consid- 
ered. 

There are several features of this plot that must be not- 
ed. The most important aspect is that the curves for ice 
thicknesses of 0.5, 0.25, and 0.125 inches are identical. 

This shows that, for these thicknesses, the outer convective 
surface has no influence on the time it takes for the ice- 
abrasion shield interface to reach 32° F. The curve corre- 
sponding to this set of ice thicknesses will be called the 
"base curve." 
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Two other important features will be called the "critical 
thickness" and the "deviation temperature rise." The criti- 
cal thickness is the minimum thickness of ice that will not 
result in a deviation from the base curve. In Figure 10, 
the critical thickness is approximately 0.125 in. For ice 
thicknesses above 0.125 inches, the temperature rise versus 
deicing time curve will not deviate from the base curve. 
For ice thicknesses below 0.125 inches, the curve will devi- 
ate from the base curve at the deviation temperature rise. 
The deviation temperature rise for 0.0625 inches is 42 °F, 
and for 0.0313 inch thick ice it is 37°F. 

Two points from Figure 10 were chosen as a way of demon- 
strating the conditions under which a variable ice thickness 
model is necessary to provide accurate results. The first 
is a point on the base curve at a thickness of 0.125 inches 
and a temperature rise of 60 °F (initial temperature at -28° 
F) . The second point is on the 0.0625 inch curve at a temp- 
erature rise of 60° F. It is important to note that the sec- 
ond point is below the critical thickness and above the de- 
viation temperature rise. It does not lie on the base 
curve . 

Figures 11 and 12 show the two-dimensional variable 
thickness ice layer results for two mean ice layer thick- 
nesses, 0.125 and 0.0625 inches, respectively. The initial 
temperature for each case is -28 °F, so that these figures 
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correspond to the two points mentioned above. The ice 
length m each case is 0.25 and 0.125 inches, respectively, 
and there is no heater gap. For each mean thickness, five 
cases of a sloped ice layer are presented, where the slope 
of the ice shape is designated by 0 , 10 , 20 , 30 , and 40.° 
The location of the melt line is indicated as a function of 
time . 

These figures clearly show the relationship between ice 
geometry, critical thickness, and deicer performance. For a 
set of conditions corresponding to a point on the base 
curve. Figure 11, where the mean ice thickness is greater or 
equal to the critical thickness, the ice geometry has a lim- 
ited effect on deicer performance. Note that even with 
drastic changes in ice shape, where a significant portion of 
the ice is less than the critical thickness, there is rela- 
tively little effect on deicer performace as is shown by the 
gradual slope of the melt curves and the relatively constant 
melt times. For these cases a one-dimensional solution 
should provide sufficiently accurate results. 

The effect on deicer performance is very pronounced when 
the conditions do not correspond to a point on the base 
curve. Figure 12, where the mean ice thickness is less than 
the critical thickness. The variation from the 0 ° case to 
the 40 0 case takes one from good deicer performance to al- 
most no deicing at all. This plainly illustrates the need 
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for a two-dimensional code capable of handling a variable 
ice layer thickness to accurately predict the thermal tran- 
sients for cases not corresponding to points on the base 
curve . 

The base curve, critical thickness, and deviation temper- 
ature rise appear to be excellent parameters for comparing 
various deicer designs. For similar power inputs and exter- 
nal conditions, these parameters will be functions of deicer 
configuration and design only. The magnitudes of these val- 
ues can be easily found using Marano's [3] one-dimensional 
code if the heater gap is not a consideration, or Chao's [4] 
two-dimensional code if the heater gap is a consideration. 
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Chapter 4 

CONCLUSIONS AND RECOMMENDATIONS 

A variable ice thickness algorithm was successfully ap- 
plied to Chao's [4] two-dimensional transient numerical mod- 
el of electrothermal deicer performance. The new model has 
been shown to be equivalent to Marano's [3] and Chao's [4] 
one and two-dimensional models for uniform ice thickness. 
It generates consistent and uniform results for variable 
thickness ice layers and allows location of the ice-liquid 
interface at any point in time. 

Determining the conditions under which a variable thick- 
ness ice layer has an effect on deicer performance proved to 
be difficult. An analysis using a one-dimensional model has 
been shown to be useful in determining these conditions. 
This analysis has shown that a one-dimensional model is suf- 
ficient for problems where the mean ice thickness is equal 
to or greater than the critical thickness. It has also 
shown that a two-dimensional analysis is necessary in cases 
where the mean ice thickness is less than the critical 
thickness. The results illustrate that the thinner the mean 
thickness of the ice, the more important the ice shape be- 
comes in order to predict the movement of the melt inter- 
face . 
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Further research in this area should be directed toward 


developing a model that can approximate the actual two-di- 
mensional blade geometry in addition to the variable ice 
thickness. A numerical model being developed at The Univer- 
sity of Toledo for this purpose is nearing completion. It 
uses a coordinate transformation technique to model the 
blade geometry, and the enthalpy method is used to model the 
phase change. 
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PART III 


EXPERIMENTAL TESTING AND RESULTS 
Chapter 5 

EXPERIMENTAL HISTORY 

Until recently, there has been very little experimental 
information available concerning the thermal behavior of an 
electrothermal deicer. Commercial manufacturers often poss- 
ess such information, but consider it proprietary and do not 
publish it in the open literature. The existing available 
information comes from both laboratory and flight test data. 

Gent and Cansdale [10] have published experimental tran- 
sient temperature profiles from laboratory deicer pads. 
They used these results to validate their one-dimensional 
computer model. Marano [3] also used these results for the 
same purpose. The signifigance of this data is limited be- 
cause the deicer pads tested were flat plate models of full 
size helicopter blade-deicer constructions. Actual deicers 
may not respond in the same manner during flight conditions. 
Flight test or wind tunnel data is thought to be more perti- 
nent . 

An early attempt at obtaining flight test data involved a 
UH-1H helicopter equipped with an electrothermal deicing 
system on the main and tail rotors [17]. Blade surface 
temperature was recorded during dry and icing conditions. 
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Transient data was not established because only peak temper- 
atures were recorded during the deicing cycle. This infor- 
mation was not detailed enough to allow a rigorous evalua- 
tion of one and two-dimensional numerical models. 

A very recent and detailed investigation promises to pro- 
vide valuable flight data. Electrothermal deicing tests of 
a Royal Air Force HC-Mkl Chinook helicopter have been con- 
ducted (18] by the Boeing-Vertol Company. One of the blades 
was equipped with both internal and external temperature 
sensors. Tests were run during both dry and icing condi- 
tions. A preliminary analysis of the flight data using the 
Baliga [1] one-dimensional computer model was conducted. 
The Baliga model tended to overpredict the actual tempera- 
tures, but the overall correlation appeared to be quite 
good. 

Due to this lack of avilable experimental data which was 
needed to evaluate The University of Toledo computer codes, 
the tests described in the following chapter were conducted. 
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Chapter 6 

FORMULATION OF THE EXPERIMENTAL METHOD 


6.1 OBJECTIVE 

During the past few years, icing research sponsored by 
the NASA Lewis Research Center has yielded several analyt- 
ical computer programs developed for the prediction of the 
thermal behavior of an electrothermal deicer. Both one and 
two dimensional models have been developed, and different 
approaches have been used. It has not been possible, how- 
ever, to validate these codes because an experimental data- 
base for electrothermal deicer performance has never been 
established. Therefore, a main objective of the experimen- 
tal testing was to develop an experimental database, there- 
by providing a standard against which the computer codes 
could be compared. In so doing, the limitations of the ma- 
thematical models could be determined. 

Another objective of this testing was to further document 
the ice accretion process and the resulting ice shapes. A 
computer program is being developed at NASA Lewis that will 
allow the prediction of an ice shape for a given set of con- 
ditions. Information collected during the testing would be 
useful m validating this numerical model of the ice accre- 
tion process. 
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Finally, a more general directive of the test was to re- 
cord the thermal response of the blade during two different 
sets of conditions. The first was a condition of evapora- 
tive cooling while the blade was running wet at temperatures 
above the freezing point. The second was an icing condition 
at temperatures below the freezing point. The entire ther- 
mal response, beginning with the onset of the spray, was to 
be recorded. This information would enable quantitative in- 
terpretations to be made concerning the physics occurring 
during the icing and deicing operations. 

Accomplishing these objectives would go far in establish- 
ing an experimental database documenting the thermal behav- 
ior of an electrothermal deicer. It would also provide much 
information concerning the thermodynamic response of an air- 
foil to icing and convective cooling conditions. 

6 . 2 TEST PLAN 

A UH-1H (NACA 0012) helicopter blade was outfitted with 
an electrothermal deicer designed and manufactured by the B. 
F. Goodrich Company. The resulting composite blade was 
heavily instrumented with thermocouples during the lamina- 
tion process. The model was then mounted on the oscillating 
rig in the test section of the Icing Research Tunnel (IRT) 
at NASA Lewis. Data acquisiton consisted of subjecting the 
blade, which was held stationary in the present tests, to a 
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wide range of aerodynamic and icing conditions and recording 
the response. 

The test conditions necessary to satisfy the objectives 
resulted in a four phase test program. The first phase was 
a dry test with no spray. Runs both with and without heat 
were recorded. These runs were meant to provide information 
on thermocouple consistency and performance, and to document 
the blade thermal behavior under tunnel flow conditions. 

The second phase consisted of running the blade wet at 
temperatures above freezing, both with and without heat. 
The entire process temperatures, from the onset of the 
spray, were recorded for all runs. Those runs with heat 
were meant to provide information on the thermal response of 
the deicer while running wet. Those tests without heat were 
meant to record the effect of convective cooling on the 
blade temperature. 

Phase three was an ice accretion documentation test. The 
heaters were not activated, which allowed the operators to 
enter the tunnel and record the ice shape after the ice ac- 
cretion process was completed. Both tracings and photo- 
graphs of the ice shapes and their cross sections were re- 
corded. This phase was vital because the ice shapes would 
be lost during the deicing tests. 
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The final phase comprised the deicing tests. Icing con- 
ditions during these tests were meant to match those during 
the ice accretion documentation tests. This would allow po- 
sitive identification of the ice shapes lost during shed- 
ding. 

There were two parts to the deicing tests. The first 
part involved a simultaneous cycling of the eight heater 
zones, in which all of the heater zones were turned on and 
off at the same time. The second part consisted of a phased 
cycling of the heater zones. Each zone was cycled with dif- 
ferent turn-on and turn-off times. The data analyzed in 
the present study was taken during the simultaneous cycling 
part of the testing. All of the testing was done during the 
period from January 31 to March 1, 1985. 


6.3 INSTRUMENTATION 

The model construction started with a section of a UH-1H 
(NACA 0012) helicopter blade. The blade was mostly aluminum 
with a stainless steel abrasion shield, as shown in Figure 
13. The electrothermal deicer consisted of a layer of epox- 
y/glass insulation, a copper heater element, another layer 
of epoxy/glass insulation, and a stainless steel abrasion 
shield. The layers of the deicer were laminated to the 
blade using epoxy adhesive. A detailed description of the 
UH-1H blade and electrothermal deicer construction is pre- 
sented in Figure 14 and Table 2. 
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The heater consisted of eight individual zones, 1.0 inch 
wide, which were electrically insulated from each other. 
Each zone was divided into a winding sub-element, 0.165 
inches wide, with a gap 0.061 inches wide, which was neces- 
sary for electrical insulation. Subdividing the zones re- 
sulted in a more uniform heater output. A single heater 
zone layout is illustrated in Figure 15. The heater element 
was manufactured by B. F. Goodrich using an etching process, 
and was attached to a backing material that prevented dam- 
age . 

Instrumentation consisted of afixmg thermocouples to the 
blade and deicer in three layers and at three planes, as 
shown in Figures 16, 17, and 18. The planes were labeled 
hh, BB, and CC, corresponding to an upper plane, a middle 
plane, and a lower plane, respectively. The three layers of 
thermocouples were placed at the inner side of the D-spar, 
the inner side of the heater mat, and the outer side of the 
abrasion shield. 

The inner side of the D-spar had eight thermocouples po- 
sitioned one inch apart as closely as possible to the center 
of the heater zones. Thermocouples on the mat were posi- 
tioned one inch apart in the center of each of the eight 
zones, and were attached to the center of the heater subele- 
ments rather than over a gap, as shown m Figure 15. A to- 
tal of eleven thermocouples were attached, one inch apart, 


68 



to the outer surface of the abrasion shield. Eight were 
centered over heater zones, and the remaining three were 
spaced over an unheated portion of the blade. 

All of the thermocouples were constructed of chromel and 
alumel wire, 0.003 inches m diameter. Those at the heater 
and substrate were spot welded at the junction and then at- 
tached to the desired location. The abrasion shield thermo- 
couples were fabricated by drilling a hole through the abra- 
sion shield, inserting the thermocouple wires, spot welding 
the junction, and dressing the surface. 

After fabrication, the thermocouples were calibrated at 
two temperatures. The first was at room temperature, 74 °F, 
and the second was at 32 °F using an ice bath. In both cas- 
es, the thermocouples read within ±1°F. 

After a few preliminary bench tests, it became apparent 
that many of the thermocouples were inaccurate. The thermo- 
couples were very fragile, and several had been damaged dur- 
ing the lamination process. Figures 16, 17, and 18 indicate 
the twenty-seven thermocouples that were disregarded as a 
result of this initial testing. These thermocouples were 
not used later in this study. 

During the testing, the thermocouples were interfaced 
with the Escort II data acquisition system. This system 
considerably simplified the data acquisition, allowing large 
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numbers of thermocouple readings to be recorded quickly and 
efficiently. The Escort II system was operated in superplex 
mode, and initially scanned the fifty-four thermocouple 
readings once every 0.3 seconds. After operating the system 
for a few initial tests, it became apparent that 0.3 seconds 
per scan was producing more data than could be efficiently 
processed. The scan time was increased to 0.9 seconds, aft- 
erwhich no further problems were encountered. 

A special computer program, IRT1D009, was developed that 
transferred the readings from Escort to an IBM 370 computer. 
Once the readings were transferred to the 370, they were 
copied onto tape for permanent storage. The computer pro- 
gram also processed several different types of graphical 
plots of the experimental data. This program was capable of 
storing, retrieving, and modifying the experimental data- 
base . 

The heater zones for the simultaneous cycling portion of 
the deicing tests were connected to three different power 
supplies. The top three zones, the zone below the stagna- 
tion point, and the bottom four zones formed the three inde- 
pendent groups. This arrangement allowed different power 
densities to be applied to different portions of the heater. 
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6.4 DETERMINING THE TEST CONDITIONS 


As mentioned previously, the testing was divided into 
four phases. The test matrix for each phase is presented m 
Tables 3 through 6. The definitions of the symbols repre- 
senting the conditions are as follows: 

= test section air velocity (miles/hr.); 

T = tunnel air stagnation temperature (°F); 

00 

a = blade angle of attack (degrees); 

PI = power supplied to bottom four heater zones 
(Watts/m. 2 ) ; 

P2 = power supplied to top three heater zones (Watts/ 

m. 2 ) ; 

P3 = power supplied to heater zone just below the 

2 

stagnation point (Watts/m. ); 
t = heater on time (sec.); 
t = heater off time (sec.); 
n = number of heater cycles; 

d = droplet diameter (microns); 

3 

LWC = liquid water content (g/m ); and 

t . = time of ice accretion (min.), 

ice 


Selection of the test condition matrix was somewhat arbi- 
trary. Duplication of a set of test conditions was avoided. 
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and the individual conditions were varied as much as possi- 
ble. The most stringent criterion was the matching of icing 
conditions between the ice accretion documentation tests and 
the deicing tests. 

6 . 5 TEST PROCEDURE 

6.5.1 Phase I - Dry Air Tests 

The procedure for the dry air tests was: 

1. Establish the desired conditions for V , a , and T ; 

2. Allow the blade to reach equilibrium at T^ ; 

3. Set the on/off times and power densities for the 
eight heater zones; 

4. Activate the Escort data acquisition system and re- 
cord the initial temperature distribution; 

5. Turn on the heaters; and 

6. Turn the heaters off, deactivate Escort, and return 
the tunnel to idle speed. 

Steps 3 and 5 were skipped if an unheated run was de- 
sired. 

Initially, the heaters were turned off if they reached 
temperatures greater than 160 6 F . After a few tests, it was 
found that they could tolerate much higher transient temper- 
atures. An upper limit of about 200°F was then allowed be- 
fore cutting power to the heater zones. There were a total 
of 46 dry runs, 35 being with heat applied, and 11 without 
heat . 
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6.5.2 Phase II - Wet Air Tests 


The test procedure for the second phase was: 

1. Establish the desired conditions for V , a, and T ; 

2. Allow the blade to reach equilibrium at T^ ; 

3. Set the on/off times and power densities for the 
eight heater zones; 

4. Activate the Escort data acquisition system and re- 
cord the initial temperature distribution; 

5. Turn on the cloud spray with the desired conditions 
of LWC and d; 

6. Record the temperature response of the blade for five 
minutes ; 

7. Turn on the heaters; and 

8. Turn the heaters off after the desired number of cy- 
cles, deactivate Escort, and turn the spray off. 

Steps 3 and 7 were skipped if an unheated run was desired. 
In this case. Escort was deactivated after equilibrium was 
reached. There were a total of 45 wet runs, 10 runs being 
with heat applied and 35 runs being without heat. The same 
temperature limit for the heater zones was enforced as for 
the Phase I tests. 

6.5.3 Phase III - Accretion Documentation Tests 

The ice accretion documentation tests were performed as 
follows : 
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1. Establish the desired conditions for V OT , a, and T^ ; 

2. Allow the blade to reach equilibrium at ; 

3. Activate the Escort data acquisition system and re- 
cord the initial temperature distribution; 

4. Turn on the cloud spray with the desired conditions 
of LWC and d; 

5. Allow the ice accretion process to continue for five 
minutes from spray activation; 

6. Turn off the spray and deactivate Escort; 

7. Return the tunnel to idle speed; 

8. Enter the tunnel and document ice accretion by trac- 
ing profiles and photographing ice cross-sections; 
and 

9. Remove the remaining accreted ice. 

The accreted ice profiles were traced by sectioning the 
ice shape with a steam knife and inserting a cardboard temp- 
late into the slot. A line was then drawn around the shape 
onto the cardboard. This was sufficient for recording the 
general ice shape, but fell short with regard to the details 
of the ice structure. 

Another method was developed because of the presence of 
the electrothermal deicer. The deicer was turned on just 
long enough to destroy the ice adhesion at the abrasion 
shield surface. The ice shape was then carefully removed and 
sectioned with a steam knife. These sections were then pho- 
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tographed next to a scale m order to provide a dimensional 
reference. Almost all of the ice structure detail was pre- 
served if enough care was exercised in its removal. This 
process resulted m some interesting and unique photographs. 
There were a total of 30 ice accretion documentation runs. 

6.5.4 Phase IV - Deicing Tests 

The procedure for the simultaneous heater zone activation 
tests was: 

1. Establish the desired conditions for V , a , T • 

00 00 

2. Allow the blade to reach equilibrium at T^ ; 

3. Set the on/off times and power densities for the 
eight heater zones; 

4. Activate the Escort data acquisition system and re- 
cord the initial temperature distribution; 

5. Turn on the cloud spray with the desired conditions 
of LWC and d; 

6. Record the temperature response of the blade for t ^ ce 
minutes; 

7. Turn on the heaters; and 

8. Turn the heaters off after the desired number of cy- 
cles, deactivate Escort, and turn the spray off. 

The same heater temperature restrictions as for the Phase I 
and II tests were enforced in order to prevent heater burn- 
out . 
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An identical procedure was used for the phased heater 


zone activation tests except for step 3, which was; 

3. Set zonal phase lags, on/off times, and power densi- 
ties. 

There were 72 simultaneous heater cycling tests and 27 
phased heater cycling tests. 
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Chapter 7 

DISCUSSION OF RESULTS 


7.1 PHYSICAL DESCRIPTION AND ANALYSIS OF THE EXPERIMENTAL 
DATA 

A vast amount of data was acquired during the testing. 
There were a total of 211 tests over a wide range of condi- 
tions. A breakdown as to the types and numbers of the tests 
is included in the test matrix. For obvious reasons, it was 
not possible to analyze every one of the tests. It was nec- 
essary to reduce the data to an analyzable set. This was 
accomplished in two steps. A preliminary set of about 30 
readings was randomly selected from the database. This set 
was studied in order to get a feel for the general physical 
trends that were occurring in the data. The second step was 
a further reduction to ten readings, allowing a detailed 
comparison to the numerical simulation models. The tempera- 
ture versus time curves from these ten readings for selected 
thermocouples are shown in Figures 19 through 28. Figures 
19 through 23 represent the dry heat tests, and Figures 24 
through through 28 are the deicing tests. Many of the gen- 
eral physical characteristics and trends about to be dis- 
cussed are shown on these figures. 
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7-1.1 General Physical Description 
7. 1.1.1 Dry Runs 

The abrasion shield temperature response for the dry runs 
was observed to be very dependent upon the external flow 
conditions. This was expected since the external heat 
transfer coefficient was relatively large at the outer sur- 
face. Because of this, the abrasion shield was always at 
the flow freestream stagnation temperature before heat was 
applied. 

In contrast, the initial readings of the heater and D- 
spar thermocouples were never at the flow freestream stagna- 
tion temperature. This behavior resulted from not allowing 
the blade to cold soak for long enough periods of time be- 
tween runs. Residual heat from previous runs had a tendency 
to remain in the blade. A low convection coefficient at the 
inner surface coupled with the relatively large mass of the 
D-spar caused this condition. 

In general, the temperature response at any layer was 
strongly influenced by the heater power density and on/off 
times. Most of the influence was in maximum temperature 
magnitude. At a given layer, the shapes of the temperature 
versus time curves were quite similar for different posi- 
tions and conditions. 
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7. 1.1.2 Deicing Runs 


The thermal response at the abrasion shield for the deic- 
ing runs was visibly different from that for the dry runs. 
The initial temperature was seldom at the freestream stagna- 
tion temperature or consistent from cross-section to cross- 
section. This behavior was caused by a fluctuation in temp- 
erature at the onset of ice accretion. 

Two types of fluctuations occured. A rise in temperature 
due to a latent heat effect, or a dip in temperature due to 
convective cooling. A correlation relating the type of 
fluctuation to icing conditions or position was not found. 
In either case, the largest deflection occurred at the onset 
of ice accretion, afterwhich the temperature tended to ap- 
proach the ambient. 

The D-spar and heater behaved in much the same manner as 
for the dry runs. At the heater, some of the readings 
showed a slight temperature deflection at the onset of ice 
accretion. The effect was similar to, but much more limited 
than at the abrasion ahield. The D-spar never showed this 
behavior and appeared to be quite independent of the outer 
surface conditions. 
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7.1.2 Reducing the Database to an Analyzable Set 

Narrowing the database down to an analyzable set that 
could be compared to the numerical models was necessary. A 
total of ten runs was chosen for the final database. This 
number was completely arbitrary. Five dry runs and five 
deicing runs were chosen. A set of criteria was used m se- 
lecting the ten readings. The deicing runs were chosen 
first using the following guidelines: 

1. The runs were limited to a 0° angle of attack; 

2. As wide a variety of ice shapes as possible was de- 
sired; 

3. Runs chosen were limited to those for which the ice 
shape was definitely known; and 

4. If possible, runs where shedding occurred were cho- 
sen. 

Based on the set of deicing runs chosen, a set of dry 
runs with heat was chosen using the following restrictions: 

1. Limit to a 0° angle of attack; 

2. Attempt to match the conditions as closely as possi- 
ble with those of the deicing runs already chosen; 

3. Try to -get a wide range of temperatures; and 

4. As a last resort, vary the heater power densities. 

The ten readings chosen and the conditions for each are 
shown in Table 7. Figures 19 through 28 present the temper- 
ature profiles for these cases. A general overview of the 
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temperature profiles of this set showed that it was repre- 
sentative and consistent with the physical trends already 
noted. 

7.1.3 Thermocouple Response as Related to Consistency 

An important part of validating the thermocouple readings 
was to check for consistency. Behavior should be similar 
from reading to reading or position to position. If a ther- 
mocouple exhibited erratic behavior for any of the ten read- 
ings, it was discarded. 

For example, thermocouples 59 and 64 exhibited consistent 
behavior for the dry runs but not for the deicing runs. 
Thermocouple 59 did not respond at all in the deicing runs, 
and thermocouple 64 did not appear to behave consistently 
with respect to those at other positions. This was consid- 
ered erratic behavior. Thermocouple 75 did not respond at 
all m the deicing runs, and thermocouple 76 did not respond 
in the dry runs and appeared to give inconsistent readings 
in the deicing runs. 

The judgement as to whether or not a thermocouple was 
consistent or erratic was somewhat arbitrary. Disregarding 
a thermocouple was not usually critical because of the pres- 
ence of duplicates for each position. It was safer to just 
disregard a thermocouple if there was any doubt about its 
validity. This had an effect only on the deicing tests, and 
involved thermocouples 64 and 76, which were disregarded. 
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7.1.4 Thermocouple Response as Related to Cross-Section 
Dependency 

7. 1.4.1 Dry Runs 

To aid m analyzing the data, "cross-section dependency 
between planes AA, BB, and CC was tested by plotting all 
three thermocouples for a given position on one graph. Some 
of the positions had only one good thermocouple. 

As was mentioned previously, the calibration showed that 
the thermocouples were accurate within ± 1°F. This made the 
experimental error range 2°F for the dry runs. Thermocouple 
readings were considered to be cross-section independent if 
they were within 2 °F of each other at a given position for 
all three cross-sections. 

Most of the thermocouple readings at the abrasion shield 
and the D-spar were within 2°F of each other for the three 
cross-sections at a given position. This behavior indicated 
that the thermocouple readings were cross-section indepen- 
dent for the abrasion shield and D-spar. There were no du- 
plicates at the heater, so this analysis did not apply. 

7. 1.4.2 Deicing Runs 

At first, the deicing runs appeared to be cross section- 
dependent at the abrasion shield. The peak magnitude of the 
thermal response varied from cross-section to cross-section 
at a given position. After further study, however, the 


82 



cross-section discrepancies appeared to be caused by differ- 
ences in initial temperature. These differences were caused 
by the onset of ice accretion, as discussed m a previous 
section. A different approach was necessary to analyze the 
deicing runs. 

The new approach used the temperature difference from 
initial to peak temperature m the cross-section compari- 
sons. A temperature rise parameter, AT, was defined as the 
total rise from initial to peak temperature. Using this pa- 
rameter m the cross-section comparisons removed the initial 
temperature dependency. 

Before AT could be compared amongst the cross-sections, 
the experimental temperature error range needed to be estab- 
lished. As for the dry runs, the accuracy of any tempera- 
ture measurement was within ±1°F, so that the error range of 
each measurement was 2°F. Therefore, the experimental temp- 
erature difference error range was 4°F. 

Most of the temperature differences at the abrasion 
shield for the three cross-sections were within 4° F of each 
other at a given position. This indicated that the abrasion 
shield thermocouples were cross-section independent. In 
other words, the response at one cross-section would be rep- 
resentative of the other two cross-sections for a given po- 
sition. 
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This analysis was not necessary for the heater and the 
D-spar. There were no duplicates at the heater and the ice 
accretion temperature fluctuation did not penetrate to the 
D-spar. The D-spar behaved exactly the same as in the dry 
runs, with a few exceptions. 

These exceptions were due to the presence of one inch 
wide anti-icers at the top and bottom of the heater mat. 
Their presence was meant to prevent ice bridging, which 
could delay or stop ice from shedding. After long periods 
of time, the initial temperature at cross-sections AA and CC 
began to rise above cross-section BB. Heat from the anti-i- 
cers was being conducted to the D-spar. Eventually, it 
reached the thermocouples at cross-sections AA and CC. This 
phenomenon was intensified by the high thermal conducitivity 
of the aluminum. 

Overall, a minority of readings showed this effect. Most 
of the substrate thermocouples behaved very similar to the 
dry runs. 

7. 1.4.3 Establishing a Representative Set of Thermocouples 

Having shown that the thermocouple readings were cross- 
section independent allowed the selection of a representa- 
tive group. This was an important step as it eliminated 
considering cross-sections during the numerical model vali- 
dation. 
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The representative set of thermocouples is shown in Table 
8. The middle cross-section, BB, was chosen at the D-spar. 
This cross-section was least affected by heat conducted from 
the anti-icers for the deicing runs. There was no choice at 
the heater because there were no duplicates. At the abra- 
sion shield, the top cross-section was chosen if possible. 
The middle cross-section was chosen if there was not a good 
thermocouple at the top cross-section for a given position. 
These decisions were not extremely critical since the ther- 
mocouple response had been shown to be cross-section inde- 
pendent . 

7.1.5 Thermocouple Response as Related to Position 

An analysis studying the relationship between thermocou- 
ple response and position was in order. The object of this 
analysis was to look for and explain any physical trends m 
relation to position that were occurring in the experimental 
data . 


7 . 1 . 5 . 1 Dry Runs 

A definite positional trend was noted at the abrasion 
shield for the dry runs. It was noted to occur for all of 
the readings. A technique was used to bring this trend onto 
the same order of magnitude for all five of the dry read- 
ings . 
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The temperature rise, AT, from initial to peak tempera- 
ture was measured and tabulated for each position and every 
run. They were brought onto the same order of magnitude by 
dividing by the average temperature rise, (AT)avg, for each 
run. This technique produced surprisingly consistent re- 
sults, as may be seen m Figure 29. 

This plot shows a physical trend consistent for all of 
the dry runs. Positions 1 through 3 show an increasing mag- 
nitude in AT. There are no heaters at these positions, and 
lateral heat conduction is solely responsible for the in- 
creasing trend. Position 4 shows a sudden increase in AT 
due to the start of the heater zones. Positions 4 through 
11 have heater zones. 

Starting at position 4, the curve slopes downward to a 
minimum at position 6, the stagnation point. The curve then 
starts back up and increases from positions 6 through 8. 
The curve is almost symmetrical about a cusp at the stagna- 
tion point. This was caused by the flow symmetry and a con- 
vection coefficient which increased to a maximum at the 
stagnation point. The UH-1H blade is a symmetrical airfoil 
and there was a 0° angle of attack. 

Progressing from position 8 to 9 shows a dip m the 
curve. This was possibly caused by an air gap running the 
length of the blade somewhere between the heater and the ab- 
rasion shield. An air gap would impede heat flow to the ab- 
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rasion shield and result in a lower peak temperature. How- 
ever, it is also possible that this dip could have been a 
consequence of the flow behavior around the blade. 

The curve continues to rise from position 9 to 11. This 
corresponds to the further development of the boundary lay- 
er. Correspondingly, for a flat plate in parallel flow, the 
convection coefficient decreases with increasing distance 
from the leading edge. This trend occurs whether the flow 
is laminar or turbulent. 

The curve in Figure 29 is a physical characteristic of 
the blade. This hypothesis is supported by the agreement of 
the curve for all five of the dry runs. Heater on-times, 
power densities, and external flow conditions change, but 
the curve remains similar from reading to reading. 

The same analysis used for the abrasion shield was car- 
ried out for the heater thermocouples. Once again, a con- 
sistent trend from reading to reading surfaced in the exper- 
imental data, as shown in Figure 30. 

This curve shows that the temperature rise, AT, oscil- 
lates m a sinusoidal manner about the average temperure 
rise, (AT)avg, in relation to position. This behavior was 
probably caused by a combination of variations in heater 
zone power output and varying contact resistances. The con- 
tact resistances were most likely induced by air gaps in the 
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construction. The data appear to indicate that this curve 
is a physical characteristic of the blade. It changes very 
little despite changes in heater power density and on/off 
times . 

In contrast to the abrasion shield and heater, the D-spar 
thermocouple temperature profiles for each test were almost 
identical despite position. For example, a thermocouple at 
the inside of the D-spar at position 6 produced the same 
profile as a thermocouple at position 10. Even the brass 
noseblock at the stagnation point had no visible effect. 

This position independence was caused by the construction 
and physical characteristics of the D-spar. The D-spar is a 
single piece of aluminum extrusion, and aluminum has a very 
high thermal conductivity. Energy conducted to the spar was 
evenly distributed throughout almost instantaneously. At 
any point in time, a sort of equilibrium existed around the 
D-spar resulting in position independent temperature pro- 
files . 

7. 1.5.2 Deicing Runs 

The dry run analysis did not apply to the deicing runs. 
Conditions between runs varied because of the different ice 
shapes at the outer surface. A consistent physical trend 
from run to run did not occur. As a result, a different ap- 
proach was used in the analysis of the deicing runs. 


88 



The temperature rise, AT, was plotted versus position for 
each cycle. The melting point was included on the curves. 
The melting point temperature rise was calculated by taking 
the difference between the melting point, 32 °F, and the ini- 
tial temperature. Graphs for each of the deicing runs are 
presented on Figures 31 through 35. 

Figure 31 shows a set of curves from reading 197, for 
which the IRT temperature was 31 °F. Recorded absolute temp- 
eratures for this test are given on Figure 24. The tempera- 
ture rise begins at position 3. There were no heater zones 
at positions 1, 2 or 3 , so this temperature increase at po- 
sition 3 was entirely due to lateral heat conduction. 

Positions 5, 6, and 7 are in the water droplet impinge- 
ment zone. Prior to heater activation, there was a glaze 
ice shape over these regions. The ice cap was shed before 
the peak temperature was reached. It did not reform before 
the next cycle because the abrasion shield was too warm. 
Position 6, the stagnation point, was not the point of low- 
est temperature as it was in the dry runs. Figures 24e, f 
and g illustrate the melting and shedding that occurred at 
positions 5, 6 and 7 at 32°F. The slight change in slope 

that is evident in the abrasion shield temperature just af- 
ter the heater is turned on is characteristic of this phe- 
nomena . 
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Moving on to positions 8, 9, and 10 in Figure 31 shows a 
steady increase in temperature rise. This was probably due 
to the development of a turbulent boundary layer and de- 
creasing convection coefficient. 

Reading 209 is shown on Figure 32, with the absolute 
temperatures given on Figure 25. The IRT temperature for 
this reading was 28°F. This test is similar to reading 197. 
An ice cap was formed before each of the three cycles. How- 
ever, it was shed before the temperature peak at the abra- 
sion shield was attained for all three cycles. The same 
trends and physical explanations apply as in reading 197. 
Positions 1 and 2 show evidence of lateral conduction. Fig- 
ure 25f at position 6 clearly illustrates that melting and 
shedding occurred at 32 °F during the first two cycles. Fig- 
ure 25e at position 5 illustrates that melting and shedding 
occurred during all three cycles. The response at position 
7, Figure 25g, does not show a clear enough change in slope 
to indicate these effects. 

Reading 213 is also similar to reading 197, and is pre- 
sented on Figure 33. The IRT temperature for this reading 
was 24 °F. Absolute temperatures are given on Figure 26. 
Ice was formed prior to each of the three cycles but was 
shed before the temperature peaked. The ice covered posi- 
tions 6 and 7 on the blade. The physical explanation of the 
trends is the same as in reading 197. Figures 26f and g 
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show that melting and shedding occurred at positions 6 and 7 
at 32°F during the first cycle. 

Reading 234 is different from the previous readings in 
that shedding did not occur. The ice continued to build and 
increase in thickness during all three cycles. The first 
cycle in this test was 20/60, and the next two cycles were 
30/60. The approximate ice shape is shown m the tempera- 
ture rise versus position plot. Figure 34, and the absolute 
temperatures are given on Figure 27. The IRT temperature 
reading was 16° F. Starting at position 1 on Figure 34, the 
temperature increases to position 3. There were no heaters 
at these positions, and the rise was entirely due to lateral 
heat conduction. 

An interesting feature on this plot is the decrease in 
temperature rise from position 3 to position 4. There is a 
heater at position 4, but not at position 3. There is no 
immediate explanation for this behavior. 

Positions 5 through 7 are beneath the ice cap. The temp- 
erature rises for these positions are almost symmetrical 
about the stagnation point.' Position 6 is not a minimum for 
the heated zones because it is not a convective interface as 
it was in the dry runs. It should be noted that melting ap- 
peared to occur at positions 5 through 7 for all three cy- 
cles, and yet shedding did not take place. Moving outside 
the ice cap from position 8 to 10 shows an increasing temp- 
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erature rise. This is characteristic of a developing turbu- 
lent boundary layer and a decreasing convection coefficient. 
Figures 27e, f and g are interesting in that they illustrate 
that melting and refreezing occurred at positions 5 and 7 
during the second and third cycles, whereas position 6 does 
not show the required change in slope, even though the temp- 
erature exceeded 32 °F. 

Reading 275, Figures 35 and 28, was like reading 234 in 
that shedding did not occur. The IRT temperature was -3° F. 
Only at position 5 during the third cycle did melting and 
refreezmg occur. This is evident from Figures 35 and 28e. 
A rime ice shape was formed prior to the first heater cycle 
and continued to grow during the second and third cycles. 
Positions 4 through 8 were covered by the ice cap. 

Lateral heat conduction effects are evident from posi- 
tions 1 through 3 in Figure 35. A lack of symmetry under 
the ice cap is apparent. This might have been caused by 
nonuniform ice adhesion. Imperfect ice adhesion could cause 
a contact resistance, thereby impeding the rate of heat 
transfer to the ice. The high degree of anti-symmetry may 
have also been caused by inconsistencies in the ice physical 
properties. Rime ice is a somewhat porous substance. Dif- 
ferences in the quantity of air at different positions could 
change the thermal conductivity and heat capacity drastical- 
ly- 
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Progressing outside the ice cap from position 9 to 10 
shows an increasing temperature rise. As before, this was 
probably induced by a developing turbulent boundary layer 
and decreasing convection coefficient. 

A thin layer of frost was observed on the blade at posi- 
tions outside the water droplet impingement zone. It was 
probably condensed onto the blade by the passing moisture 
laden air. The frost appeared in a more or less random pat- 
tern, largely depending upon the ice shape. A thermocouple 
at any one position might or might not be covered by frost. 
This introduced another aspect of uncertainty when studying 
temperature profiles at the abrasion shield. 

It should be mentioned that water runback was also evi- 
dent in many of the ice accretion and deicing tests. 

Many of the dry run physical trends at the heaters were 
observed in the deicing runs. This suggested using the dry 
run analysis for the deicing runs and produced Figure 36. 
Comparing this curve to that for the dry runs. Figure 30, 
shows that they are almost identical. This behavior sup- 
ports the idea that this response is a physical characteris- 
tic of the blade. It also implies that the heater response, 
when correlated in this manner, is independent of the outer 
surface conditions. 
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The D-spar thermocouples responded in a manner almost 
identical to the dry runs. The temperature profile for each 
test was independent of position. The brass at the nose- 
block produced no discernable deviations. 

The approximate ice shapes and thicknesses for all of the 
deicing tests (197, 209, 213, 234, 275) are shown m Figure 
37. It should again be mentioned that these ice shapes were 
determined during the ice accretion documentation tests. 
The same IRT conditions were repeated, as closely as possi- 
ble, for the deicing tests. The only case where a slight 
difference occurred was in reading 197. In this case, the 
deicing test was run with a droplet diameter of 18.4 (Figure 
24), whereas the accretion test involved a droplet diameter 
of 23.1 (Figure 37) . 
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Chapter 8 

CONCLUSIONS AND RECOMMENDATIONS 

Experimental data were recorded in the NASA Lewis Icing 
Research Tunnel documenting the thermal behavior of a UH-1H 
helicopter blade fitted with an electrothermal deicer pad 
manufactured by the B. E. Goodrich Company. The testing was 
done in four phases: dry air tests; wet air tests; ice ac- 
cretion tests; and deicing tests. Originally, there were 
eighty-one thermocouples in the deicer assembly. Bench and 
consistency tests reduced this number, resulting m a final 
total of fifty-two thermocouples that were used in the anal- 
ysis. These thermocouples proved to be numerous enough and 
in such a variety of locations that an accurate record of 
the thermal response of the blade and deicer pad was ob- 
tained. 

A total of two hundred and eleven IRT readings were tak- 
en, from which ten readings, five dry runs and five deicing 
runs, were selected for further analysis. Only a zero de- 
gree angle of attack was considered. The purposes of the 
analysis were: (1) to examine the thermal response of the 
deicer assembly for information on how flow conditions af- 
fected the temperature transients (dry runs), and on how ice 
accretion and the presence of ice affected the temperature 
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profiles (deicing runs); (2) to investigate, as much as pos- 
sible, the physics of deicing; and (3) to provide a reduced 
database for comparison with mathematical models. The ther- 
mocouple readings from this set of runs were found to behave 
consistently and, within experimental error, were found to 
be independent of the three cross-sections on the blade 
where measurements were taken. This meant that the thermo- 
couple response at a given position was the same at cross- 
sections AA, BB and CC, which allowed the selection of a 
representative set of thermocouples. This further simpli- 
fied the analysis. 

For the dry runs, the magnitude of the peak temperature 
at the abrasion shield interface was found to depend upon 
position in a manner consistent for all five readings. 
Plots of AT/(AT)avg showed test independence at each blade 
position. Symmetry around the stagnation point was evident, 
as was a decreasing heat transfer coefficient with arc 
length along the blade surface. The deicing abrasion shield 
AT values were not test independent due to the different ice 
shapes, but showed moderate symmetry about the stagnation 
point. Two-dimensional effects were only noticed at the 
three blade positons that did not have heater elements. 

For both the dry and the deicing runs, the heater re- 
sponse, when correlated as AT/( AT)avg versus position, 
showed test independence for all ten runs. This behavior 
was a physical characteristic of the blade. 
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In contrast, the response of the D-spar thermocouples was 
found to be almost entirely independent of position within 
each test. Energy transferred to the D-spar was quickly 
conducted around the D-spar so that its temperature was po- 
sition independent. This behavior was caused by the high 
thermal conductivity of the aluminum, and represents a mul- 
tidimensional effect. 

In the deicing runs, melting and shedding occurred in 
three cases (tests 197 with freestream temperature 31° F, 209 
with freestream temperature 28°F, 213 with freestream temp- 
erature 24° F) before the maximum temperatures were reached. 
In reading 197, ice did not reform before the second cycle; 
for readings 209 and 213, ice reformed and shed at similar 
times during the second and third cycles as in the first cy- 
cle. For the last two deicing runs (tests 234 at 16° F, 275 
at -3 °F), the ice did not shed, although temperatures in 
reading 234 indicated that some melting had occurred. The 
occurrence of melting, refreezing , and/or shedding was docu- 
mented m the temperature response of the abrasion shield 
thermocouples . 

The tests illustrated that the criterion for shedding in 
the three cases where it did occur was that the abrasion 
shield interface temperature was 32-34°F. The ice shapes 
for these cases (glaze ice) were like frozen slush, and 
horns were evident which experienced high aerodynamic shear 
forces which aided the ice removal. 
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Test 234 at 16 °F appeared to show that melting at some 
positions occurred, but the ice shape was such that the sur- 
face stress was much less than in the previous three cases. 
The ice in this test exhibited both glaze and rime charac- 
teristics. 

As the tests were run with the primary purpose being to 
document the thermal response of the deicer pad, an evalua- 
tion of the performance of the deicer was difficult. The 

2 

power density in reading 234 was too low (8 W/m ) and the 
cycle used in reading 275 was too short (10/30), when com- 
bined with the resulting ice shapes, to obtain the necessary 
interface temperatures and conditions for shedding to occur. 

Finally, the experimental data indicated that the thermo- 
dynamic process of ice accretion had an effect on the ther- 
mal response of the deicer. This was evidenced by the in- 
flection noted to occur at the onset of ice accretion. 
However, this might have been an effect that would occur 
only in the IRT. It has not yet been determined how well 
the water spray in the tunnel emulates the supercooled water 
droplets which occur in nature. 

Further research should be directed toward studying a 
larger set of deicer runs with a wider variety of conditions 
in order to further assess deicer performance and the ther- 
mal response of the assembly. Particularly, more angles of 
attack should be considered. In addition, the oscillating 
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blade tests should be conducted since the environmental con- 


ditions would then more closely approximate in-flight condi- 
tions . 
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PART IV 


COMPARISON OF EXPERIMENTAL RESULTS WITH 
NUMERICAL CODES 

Chapter 9 

THE VALIDATION OF NUMERICAL CODES 
9.1 OBJECTIVE 

As was discussed m the Introduction, the third objective 
of the present study was to compare the experimental data 
with the numerical codes developed at The University of To- 
ledo by Baliga [1], Marano [3] and Chao [4]. This compari- 
son, or model validation, was necessary (1) in order to ver- 
ify that the model contained sufficient physics to describe 
the thermal behavior of a real deicer pad; and (2) so that 
contact resistances, etc., could be added to the model , if 
necessary, to properly characterize the construction of a 
real deicer pad, i.e., to calibrate the model for a real 
case as versus a perfectly or theoretically constructed dei- 
cer. The major result hoped for is that models will result 
that can accurately be used for electrothermal deicer de- 
sign. 

It was decided that initial comparisons of the experimen- 
tal data would be made with the one-dimensional model of 
Marano [3] in order to determine the blade locations where a 
one-dimensional model would provide accurate predictions. 
A. Peterson [18] of Boemg-Vertol has reported success with 
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the Baliga [1] one-dimensional model, but Marano's code was 
felt to be more accurate in simulating the phase change 
since it uses the Enthalpy Method to determine the location 
of the phase interface. Thus, Marano's code was chosen for 
the initial comparisons with the experimental data. 

9.2 DETERMINING THE ONE-DIMENSIONAL NUMERICAL MODEL INPUT 
PARAMETERS 

9.2.1 Modeling the Blade Geometry 

There are some obvious limitations in using a one-dimen- 
sional approximation in modeling electrothermal aircraft 
deicing. An electrothermal deicer and blade construction 
potentially constitutes a three-dimensional transient ther- 
mal problem. However, it is important to recognize which 
portions, if any, of a blade may be effectively modeled us- 
ing a one-dimensional transient numerical simulation. 

One way of looking at a one-dimensional geometry is to 
consider it as a finite thickness flat plate with infinite 
boundaries m the other two directions . If the heat gener- 
ation is uniform over the entire plate surface, energy will 
travel only inward. A boundary or discontinuity in geometry 
too close to the point of interest causes an edge effect. 
The problem is then no longer one-dimensional. A one-dimen- 
sional model is a good approximation as long as the actual 
geometry emulates a flat plate and has no end or edge ef- 
fects . 
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The flate plate analogy works best along the side of the 
D-spar, at positions 8 through 10 (Figure 19). The blade 
curvature is relatively low, and the layer thicknesses are 
constant. These positions are far enough from any disconti- 
nuities in geometry to insure minimal edge effects. There- 
fore, the best correlation should occur at this portion of 
the blade. 

The other positions depart from the one-dimensional anal- 
ogy in varying degrees. Position 11 emulates the flat 
plate geometry but loses energy to an adjacent unheated 
zone. Position 7 is right at the corner of the D-spar, and 

position 6 is affected by the brass noseblock. A somewhat 
marginal correlation should occur at these positions. 

Positions 1 through 5 were not utilized in the numerical/ 
experimental correlation. There was no heater at positions 
1 to 3, and any thermal response was due to lateral heat 
conduction only. This effect cannot be modeled using a one- 
dimensional simulation. Due to symmetry, the response for 
positions 4 and 5 was almost identical to positions 6 and 7, 
except for test 275. Modeling only positions 6 through 11 
reduced the total CPU time without sacrificing important in- 
formation. 
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9.2.2 Modeling the Composite Body 


Two cross-sections were possible depending upon the por- 
tion of the blade being modeled. Table 2a lists the physi- 
cal cross-section and properties used as input for the com- 
posite body from positions 7 through 11. At position 6, the 
brass was included, resulting in the cross-section of Table 
2b. The electrothermal deicer construction was the same for 
both cross-sections. A total of 83 nodes were used in the 
composite body without the brass, and 95 nodes were used 
with the brass. A breakdown of the number of the nodes used 
per layer is included in the tables. 

The thermal properties presented on these tables are pub- 
lished values that may vary from the actual values. Some of 
the properties were not firmly established. For example, 
the adhesive used in constructing the UH-1H blade was speci- 
fied as FM1000 film adhesive manufactured by the American 
Cyanimid Company. The thermal conductivity and heat capaci- 
ty of this substance was not known. A company representa- 
tive suggested using values for unfilled epoxy resin. 

On the other hand, many of the properties are probably 
quite accurate. Those for copper, aluminum, and stainless 
steel have been established to a high degree of accuracy. 
The only uncertainty comes from determining the grade of 
metal actually used m the blade. 
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After running a few preliminary comparisons, it became 
apparent that slightly increasing the thermal conductivities 
and heat capacities of the composite blade materials im- 
proved the numerical/experimental correlation. For all of 
the dry and deicing runs, the thermal properties input into 
the computer program were increased by a maximum of five 
percent. This was not regarded as an excessive or unrealis- 
tic approximation. 

Another parameter input into the computer simulation was 
the zonal heater power density. The heater power density 
varied from position to position due to zonal variations in 
electrical resistance, and the densities given in Tables 3, 
4, 6 and 7 are the nominal values. Actual power densities 
at each position for all ten readings are presented in Table 
9. 

9.2.3 Modeling the Ice 

The ice shapes accreted during the deicing tests were 
lost when the blade deiced. The ice thickness for a given 
position had to be known. As stated previously, this infor- 
mation was obtained by matching the icing conditions for the 
deicing tests with those from the ice accretion documenta- 
tion tests. The approximate ice shapes for the five runs 
are shown in Figure 37. The ice shapes for Readings 234 and 
275 at 2.5 minutes of accretion time were obtained by taking 
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one-half of the thickness recorded in ice accretion tests 


173 and 181, respectively, at a time of 5 minutes. This 
calculation assumed a linear rate of ice growth. 

Readings 197, 209, and 213 were ice shedding runs. Mara- 
no's [3] program handles shedding as a sudden removal of all 
of the ice and water. The ice is shed when half a specified 
node is melted. In other words, the ice is shed when the 
melt-ice interface is halfway through a particular node. 
This particular node will be called the deice node and is 
part of the program input. Its value is somewhat arbitrary. 
The thickness of a single node must be considered when se- 
lecting the deice node. 

The deice node can be any node within the ice layer. The 
abrasion shield interface is specified as node 1 of the ice 
layer. If node 1 is chosen, the ice is shed as soon as the 
interface reaches 32° F. The phase change algorithm is not 
employed. Specifying node 2 or above causes the program to 
enter the phase change algorithm before shedding. 

A large number of nodes is desirable in the ice layer. 
The enthalpy method causes a plateau at 32 °F when the ice 
melts. This occurs, as previously discussed, because each 
node is held at 32 °F until the entire nodal volume melts. 
Decreasing the nodal thickness decreases the length of time 
over which the plateau occurs. If a large number of nodes 
is used, the plateau is not descernable on a temperature 
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versus time plot. A total of ninety nodes was used m the 
ice layer. Generally, this was enough to avoid the plateau. 

There was some uncertainty m the ice thermal properties. 
The thermal conductivity of ice varies a relatively large 
amount over the experimental temperature range. Properties 
at 32° F were used because phase change occurred m most of 
the deicing runs. 

Another aspect that has an effect on ice thermal proper- 
ties is air trapped in the ice microstructure. Small bub- 
bles, a few thousands of an inch in diameter, are probably 
trapped during accretion. Some types of ice may contain 
more air than others. For all of the deicing runs, air was 
not taken into account in the ice thermal properties. 

9.2.4 Determining the Outer and Inner Surface Convection 
Coefficients 

The outer surface convection coefficients at points other 
than the stagnation point were calculated using flat plate 
correlations [19]. Both laminar and turbulent: fiow were 
considered. The equation used for laminar flow was: 


h = 0.332 - \pRe~ (Pr) 1/3 Re v < 5 x 10 5 (65) 

O X V X X 

The equation for turbulent flow was: 
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h = 0.51 (T ) 
o m 


0.3 


(1.69 pv) 

.. 0.2 


0.8 


] 


Re >2x10 


( 66 ) 


where : 


h Q = Outer surface convection coefficient, in BTU/hr- 
2 

ft -°F; 

k = Thermal conductivity of air, m BTU/hr-ft-°F; 
x = Distance along the blade surface from the stagna- 
tion point, in feet; 

Re = Reynolds number based on the distance, x, from 
the stagnation point; 

Pr = Prandtl number of the air flowing over the blade; 

= Film temperature of the air flowing over the 
blade, in degrees Rankine; 

p = Density of the air flowing over the blade, in 
lbm/f t and 

v = Local velocity of the air flowing over the blade, 
in ft/sec . 


All of the air properties were calculated at the film temp- 
erature, T m - The film temperature was approximated as the 
average between the initial and peak temperature for a given 
position at the abrasion shield. It varied from position to 
position. 
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Local velocities were used m both equations for a spe- 
cific position at the blade surface. These velocities were 
calculated using an inviscid analysis for a NACA 0012 air- 
foil [20]. The local velocities were always higher than the 
freestream velocities for the positions studied. 

The stagnation point coefficient was approximated using a 
cylinder m cross-flow correlation. The equation was taken 
from the same source [19]: 


h 


o 




0 < 6 < 60° (67) 


where : 


h Q = Stagnation point convection coefficient, m BTU/ 
hr- ft 2 - °F ; 

k = Thermal conductivity of the air flowing over the 
cylinder, in BTU/hr-f t- °F; 

D = Cylinder diameter, in feet; 

Re D = Reynolds number of the air flowing about the cyl- 
inder, based on D; and 

0 = Angle that specifies position around the cylinder 

with 0° at the stagnation point, m degrees. 
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The air properties were taken at the freestream flow stagna- 
tion temperature. The Reynolds Number was calculated using 
the freestream flow velocity, and the diameter, D, was based 
on the leading edge radius of curvature. 

Convection coefficients for the deicing runs were calcu- 
lated in the same manner. The turbulent flow flat plate 
correlation was used with air properties at the freestream 
flow temperature. Neglecting the presence of the ice shape, 
the local velocities were calculated using an inviscid anal- 
ysis for a NACA 0012 blade. The cylinder m cross-flow cor- 
relation was used at the leading edge based on the leading 
edge radius. As will be seen, these approximations produced 
acceptable results for most of the data. 

Unlike the outer blade surface, the inner surface did not 
experience forced convection. It was more like a free con- 
vective surface. Free convective surfaces usually have 

2 o 

coefficients ranging from 1 to 10 BTU/hr-ft - F. A value of 
2 

1 BTU/hr-ft -°F was used in the numerical simulation. 

9.3 COMPARING THE NUMERICAL MODEL TO THE EXPERIMENTAL DATA 
9.3.1 Dry Runs 

Having established the input parameters, the next step 
was to compare the numerical and experimental data. Mara- 
no's (3) one-dimensional model was run for all of the dry 
runs at positions 6 through 11. The results of the compari- 
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sons are shown in Figures 38 through 42. The thermocouple 
positions were modeled by plotting nodes at the spar-ambient 
interface, the inner glue-heater interface, and the abrasion 
shield- ambient interface. Included in each figure are the 
following parameters: 


= Freestream flow stagnation temperature, m °F; 

Pz = Heater power density for a specific zone, in 
W/m 2 ; 

h Q = Outer surface convection coefficient, in BTU/hr- 
ft 2 - °F ; 

t = Heater on time, in seconds; 
on 

x = Distance from the stagnation point, in inches; 
and 

Re = Reynolds number based on the distance from the 
x 

stagnation point. 


These parameters, except Re x and x , were used to generate 
the numerical data. Also included in Figures 38 through 42 
are the following symbols used to differentiate between the 
various curves: 

ABX = Experimental data at the abrasion shield; 

HX = Experimental data at the heater; 

SX = Experimental data at the D-spar; 

ABS = Numerical simulation data at the abrasion shield; 

HS = Numerical simulation data at the heater; and 
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SS = Numerical simulation data at the D-Spar. 


The Reynolds Number was an important parameter to consid- 
er. According to Reference [19], transition from laminar to 
turbulent flow starts at a Reynolds Number of 5 xlO^ and 
ends at a Reynolds Number of 2 x 10 Most of the Reynolds 
Numbers were below the transition point. This would imply 
that the flow was laminar, and that laminar coefficients 
should have been used. 

After a few trial comparisons, it become apparent that 
the laminar coefficients were too low. Likewise, the turbu- 
lent coefficients were too high. A combination of two prob- 
lems probably caused this discrepancy. First, using flat 
plate and cylindrical correlations to calculate convection 
coefficients for an airfoil is only an approximation. Sec- 
ond, tunnel turbulence probably lowered the transition Rey- 
nolds Number and raised the coefficients in general. 

Despite this discrepancy, the values used in the numeri- 
cal/experimental correlation fell between the laminar and 
turbulent values calculated using the flat plate and cylin- 
drical correlations. Convection coefficients at 65 percent 
of the fully turbulent values were used for readings 92, 93, 
and 94, and 42 percent of the turbulent values were neces- 
sary for readings 70 and 76. 
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Studying Figures 38 through 42 reveals an inconsistent 
trend in the numerical/experimental correlation. The numer- 
ical data overpredict the experimental heater temperature in 
readings 92, 93, and 94, but underpredict the experimental 
heater temperature in readings 70 and 76. This behavior can 
possibly be attributed to an error in the measured heater 
power output for readings 70 and 76. 

The input heater voltage was supposed to be at 10V for 
readings 70 and 76. At first, the power source voltmeter 
was used to set the desired voltage. However, it became ap- 
parent after a few readings that the power source voltmeter 
was not accurate enough. An external digital voltmeter, ac- 
curate to a tenth of a volt, was used to fine tune the input 
voltages for readings 92, 93, and 94, but not readings 70 
and 76. The input voltage for readings 70 and 76 was prob- 
ably higher than 10V. 

This also possibly explains why 42 percent of the fully 
turbulent convection coefficients were necessary in readings 
70 and 76 as compared to 65 percent for readings 92, 93, and 
94. A lower than actual heater power density in the numeri- 
cal model would mean less energy was reaching the abrasion 
shield. A lower external convection coefficient was neces- 
sary to match the experimental data. 

As expected, the one-dimensional model appears to be most 
accurate at positions 8, 9, and 10. At position 11, the nu- 
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merical data overpredicted the experimental data due to lat- 
eral heat conduction. Except for the heater temperature in 
Reading 94, the numerical data matched the experimental data 
surprisingly well at position 7. The effect of the change 
in geometry at the D-spar was minimized by the high conduc- 
tivity of aluminum. 

At position 6, there is a noticable breakdown m the one- 
dimensional analogy at the D-spar. The numerical data se- 
verely underpredicted the experimental temperature profile 
of the D-spar. This occurred because the brass noseblock 
included in the simulation acted as a huge heat sink. As 
discussed previously, the experimental D-spar temperature 
profiles were position independent due to lateral heat con- 
duction around the D-spar. This phenomenon could not be 
modeled using a one-dimensional simulation. 

The heater temperature was consistently overpredicted by 
the numerical data for readings 92, 93, and 94. It was 
thought that a 0.001 inch thick layer of capton, which was 
between the thermocouple and the heater, was affecting the 
measurment of the actual heater temperature. An attempt at 
modeling this geometry was made by plotting the numerical 
temperature profile at a node about 0.001 inches into the 
glue layer away from the heater. The effect on the numeri- 
cal data was found to be negligble, as shown in Figure 43. 
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9.3.2 Deicing Runs 


There were two goals to the simulation comparisons for 
the deicing runs. The first was to verify the Enthalpy 
Method as an effective model of the phase change. The sec- 
ond was to evalute the accuracy of the numerical ice shed- 
ding algorithm. 

Reading number 234 was an excellent run to use for test- 
ing the Enthalpy Method. Ice shedding did not occur, and 
the abrasion shield temperatures exceeded 32 °F under the ice 
cap. Positions 5, 6 and 7 were the only positions covered 
by ice, and positions 5 and 7 showed reasonable symmetry. 
Thus, only positions 6 and 7 were modeled. 

The composite body was modeled in the same manner as m 
the dry runs, except ice was included. The thermocouple po- 
sitions were numerically modeled by nodes at the inner am- 
bient-spar interface, the inner glue-heater interface, and 
the outer ambient-abrasion shield interface. The numerical 
and experimental comparison results for positions 6 and 7 
are presented in Figures 44 through 48 for readings 197, 
209, 213, 234, and 275, respectively. Included on these 
figures are the pertinent test conditions; 


T^ = Freestream flow stagnation temperature, in °F; 
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Pz = Heater power density for a particular position, 
2 

m W/in ; 

h Q = Outer surface convection coefficient, in BTU/hr- 
ft 2 - °F ; 

t on = Heater on time, in seconds; 

t Q ££ = Heater off time, in seconds; and 

IT = Ice thickness at a particular position . 


These parameters were important input parameters for the nu- 
merical simulation. 

Using the pertinent input parameters and the composite 
body of Table 2a with an ice layer. Figure 47, parts e 
through g, resulted for position 7 for reading 234. There 
does not appear to be any phase change in the first cycle; 
however, phase change is clearly evident in the second and 
third cycles, as was discussed above m Chapter 7, section 
7. 1.5.2. Good agreement between the experimental and numer- 
ical data resulted. The experimental results clearly illus- 
trate slope changes m the abrasion shield temperature, sig- 
nifying that melting and refreezing at 32° F did occur in the 
second and third cycles. 

2 

An outer surface convection coefficient of 145 BTU/hr-ft 
- °F was used to obtain these results. This was over twice 
that calculated using the turbulent flat plate analogy with 
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local velocities given by an NACA 0012 mviscid analysis. 
This coefficient may have been too high. The ice at posi- 
tion 7 was probably thinner than 1/16 inch, in which case a 
lower coefficient would have produced equivalent results. 

The stagnation point convection coefficient was calculat- 
ed using the cylinder in cross flow correlation with the 
airfoil leading edge radius. Using this coefficient along 
with the other pertinent input parameters produced Figure 
47, part a, for the first cycle at position 6. Even though 
the experimental data shows abrasion shield temperatures 
above 32 °F, no change in slope is evident, meaning that the 
phase change predicted by the numerical simulation did not 
occur in the experimental data. Two possible explanations 
can be offered for these results. 

The first explanation is that, considering the experimen- 
tal data to be reliable, this phenomenon was caused by an 
air gap over the abrasion shield thermocouple at position 6. 
Accreted ice is suspected as being a porous substance. 
Small air bubbles can be trapped in the ice microstructure 
during the accretion process. These bubbles are random m 
size and occurrence, causing a random effect on thermocouple 
response. Ice adhesion and density probably vary from posi- 
tion to position. 

The air effect was modeled by placing a 0.002 inch gap of 
air between the abrasion shield and the ice at position 6. 
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Including the air produced much better agreement, as may be 
seen in Figure 47, parts b through d. This figure supports 
the conclusion that air trapped in the ice during accretion 
can have a major effect on the thermal response of the abra- 
sion shield. 

The second explanation concerns the previously mentioned 
independence of sections AA, BB and CC on the blade. At po- 
sition 6, thermocouple 65 (Figure 27f) m section BB was 
chosen as the representative thermocouple since thermocouple 
54 had been eliminated by bench testing and thermocouple 76 
in section CC had been disregarded based on consistency. 
Going back to the original readings at position 6, as shown 
in Figure 49a, illustrates that if thermocouple 76 had been 
chosen, then no difference would have been observed between 
the experimental and theoretical results m Figures 47a, b 
and c. I.e., the predicted melting, which is evident in the 
response of thermocouple 76, did occur experimentally. Sim- 
ilar reasoning can be used at position 7 between thermocou- 
ples 55 (Figure 27g) and 77 (Figure 49b). Thermocouple 55 
shows melting in Figures 47, parts f and g, whereas thermo- 
couple 77 gives evidence of an air gap. 

The most plausible description of the physics that oc- 
curred m reading 234 is that section BB shows an air gap 
(thermocouple 65) at position 6 and melting (thermocouple 
66) at position 7. Correspondingly, section CC shows melt- 
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mg (thermocouple 76) at position 6 and an air gap (thermo- 
couple 77) at position 7. This would explain why no shed- 
ding occured at either cross-section in reading 234. 

As in the dry runs, the experimental D-spar temperature 
at position 6 exceeded the simulation results. 

Having demonstrated the validity of the phase change al- 
gorithm allowed the verification of the shedding algorithm. 
Readings 197, 209, and 213 were ice shedding runs. The per- 
tinent input parameters were the same, plus the following 
parameters; 

DI = Deice node; and 

MT = Melt thickness, in inches. 


As discussed previously, the deice node determines when 
shedding occurs. 

The outer surface convection coefficients were calculated 
as if the ice was not there. The flat plate turbulent cor- 
relation using local velocities from an inviscid analysis of 
an NACA 0012 was used at position 7. The stagnation point 
value, at position 6, was calculated using the cylinder m 
cross-flow correlation with the leading edge radius and 
freestream velocity. These coefficients were appropriate 
once the ice was shed. 
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Results from the numerical simulation as compared to the 
experimental data are presented in Figures 44 through 46 for 
readings 197, 209, and 213, respectively. Deicing occurred 
in the numerical data at the change in slope in the curve as 
indicated in Figure 44 for reading 197. Those simulations 
that had the deice node at node 1 did not enter the phase 
change algorithm. 

A change in slope in most of the experimental curves, as 
indicated m the plots, signifies the actual point of ice 
shedding. The test log indicated that shedding occurred at 
about ten seconds m all three runs. An accurate measure- 
ment was not taken. The numerical deicing point could be 
made to match the experimental deicing point by reducing the 
value of the deice node. 

The second cycle was modeled in reading 197. Ice did not 
accrete during the second cycle because the abrasion shield 
temperature was too high. In the second cycle, the abrasion 
shield was a convective interface and the convection coeffi- 
cient was calculated as stated above. As may be seen in 
Figure 44, parts b and d, this approximation generated a 
good correlation. 

The second and third cycles were not modeled in readings 
209 and 213 because of the complexities introduced by the 
ice accretion process. Thermodynamic effects of ice accre- 
tion were not employed in Marano's program. Besides that. 
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the ice shapes for the second and third cycles were not doc- 
umented. These discrepancies would have allowed a marginal, 
at best, correlation. 

As can be seen, the simulation results compared very fav- 
orably with the experimental results for readings 197, 209 
and 213. As in the dry runs, the experimental D-spar temp- 
erature at position 6 exceeded the simulation results in all 
cases . 

Reading 275 was the only deicing run to be numerically 
modeled all the way around the blade from position 6 through 
position 11. The outer surface convection coefficients were 
calculated using the flat plate and cylindrical correlations 
along with the NACA 0012 mviscid analysis. 

From the experimental data, no shedding occurred and 
melting is evident only at position 5 during the third cy- 
cle, as was discussed in section 7. 1.5.2 in Chapter 7. As 
may be seen in Figure 48a, the numerical simulation over- 
predicted the experimental temperature at the abrasion 
shield by a significant amount at position 6, and also pre- 
dicted a large amount of phase change in the third cycle 
which did not occur. Again, the experimental D-spar temper- 
ature was much higher than the predicted value. At position 
7, the abrasion shield prediction is much closer, but the 
heater simulation is much higher than the measured value. 
The simulation would have fit the temperatures at position 5 
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much better (Figure 28e). Positions 8 through 11 exhibit a 
good correlation, as presented in Figures 48c through 48f. 
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Chapter 10 

CONCLUSIONS AND RECOMMENDATIONS 

A one-dimensional numerical model developed by Marano [3] 
was compared to an experimental database derived from tests 
performed in the Icing Research Tunnel at the NASA Lewis Re- 
search Center. The experimental database used m the com- 
parisons consisted of thermocouple readings from five dry 
runs and five deicing runs, and was representative of a much 
larger set of readings. Input to the simulation was modi- 
fied to match experimental test conditions at a variety of 
blade positions, thereby allowing an effective comparison. 

For the dry runs, temperature responses at positions 6 
through 11 were simulated. Positions 1 through 3 had no 
heater elements, so any energy transferred to these posi- 
tions was due to two-dimensional effects. These could not 
be modeled using a one-dimensional code. Positions 4 and 5 
were also not modeled as their experimental thermal respon- 
ses showed reasonable symmetry with those from positions 7 
and 8. For the deicing runs, the temperature responses for 
the stagnation point and any positions on the lower side of 
the blade where ice had accreted were modeled. Again, be- 
cause the experimental data was taken at a zero degree angle 
of attack, the assumption of symmetry was used for the cor- 
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responding ice covered positions on the top of the blade. 
Their thermal response had shown moderate symmetry, and 
these corresponding positions were also not modeled. As was 
shown m the simulations of readings 234 and 275 m Chapter 
9, this latter assumption probably should not have been 
made, and all positions under the ice should have been simu- 
lated. 

In general, the simulation was found to overpredict the 
experimental data, especially at the heater element. This 
overprediction was possibly due to the effects of blade cur- 
vature and lateral heat conduction, neither of which were 
modeled by the one-dimensional simulation. A position at 
which the simulation consistently underpredicted the experi- 
mental data occurred at the inner side of the D-spar beneath 
the stagnation point. This behavior was caused by a break- 
down in the one-dimensional analogy at this point. 

Flat plate and cylindrical correlations were used, with 
an inviscid analysis of a NACA 0012 airfoil, in calculating 
the outer surface convection coefficients. Coefficients 
falling between the laminar and turbulent values were used 
for the dry runs, and fully turbulent values were used for 
the deicing runs. A free convection coefficient was used at 
the inner blade surface. These values generated adequate 
results m most cases. 
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The dry run simulation comparisons achieved moderate suc- 
ss. The abrasion shield temperatures showed good agree- 
nt, as did the substrate comparisons at positions other 
an the stagnation point. The predicted heater tempera- 
res were too high, as stated, but these temperatures were 
e most sensitive ones with respect to the power intensity, 
r predictive purposes, using the flat plate and cylinder 
at transfer coefficients, it appeared that the one-dimen- 
anal model can be used to model the thermal response of 
i deicer assembly. The experimental data did not appear 
show any other two-dimensional effects for positions 6 
-ough 11. It would also be expected that the model would 
jdict better at the higher ambient temperatures since the 
irmal gradients were less. Runs 70 and 76 were at 17° F 
l -4°F, and were modeled quite well except for the heater 
position 6, whereas runs 92, 93 and 94 were all for am- 
■nt temperatures of -15°F. For these latter cases, agree- 
t was not as good. No indications became evident during 
comparison process that imperfect contact between layers 
present in the deicer assembly. This would have neces- 
ated the addition of contact resistances to the simula- 
n code. 

For the deicing runs, the Enthalpy Method was shown to 
quately model the phase change which occurred at the ab- 
lon shield - ice interface. Evidence has been presented 
ch indicates that air gaps may have occurred at this m- 
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terface in reading 234, thereby delaying phase change . 
preventing ice shedding. The onset of phase change and, 
refreezing was discernable m the temperature versus t 
plots as a sudden change in slope at 32 °F. 

The Enthalpy Method and the deicing criteria set forth 
Marano's code were found to be adequate m modeling t 
first cycle of ice shedding. Experimental ice shedding c 
curred in readings 197, 209 and 213 when the abrasion shie 
- ice interface temperature was 32-34° F. Subsequent eye] 
were not modeled because an ice accretion model was not 3 
eluded m the simulation. Ice shedding was discernable 
the experimental and numerical temperature versus time pic 
as a sudden increase m slope following change m phas 
The melt thickness at which shedding occurred in the simul 
tion was always less than 0.005 inches. As in the dry rur 
no contact resistances had to be added to the numeric 
code . 

In general, the one-dimensional code of Marano show 
good comparison with the experimental data, with the compa 
isons being better at the higher freestream temperatures 
This code definitely contains sufficient physical informa- 
tion so as to adequately model the thermal response of the 
electrothermal deicer assembly at positions where two-dime 
sional effects are small. 
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Further research should be directed toward three areas. 


The first effort should be spent toward expanding the pres- 
ent investigation to include more deicing run experimental 
data constituting a wider set of conditions. The question 
of cross-section independence of the thermocouple readings 
for the deicing runs should be re-examined, and all blade 
positions that are covered by the accreted ice should be 
modeled. Two-dimensional simulations developed at The Uni- 
versity of Toledo should also be employed and evaluated. 
The second effort should be directed toward modifying Mara- 
no's code to include an internally calculated heat transfer 
coefficient at the outer surface, and also, if possible, an 
ice accretion model. Finally, the third effort should be 
set toward acquiring the flight test data compiled by the 
Boeing-Vertol Company during tests of an RAF HC-Mkl Chinook 
helicopter equipped with electrothermal deicers on the main 
rotors. This data would provide invaluable information m 
regard to the differences between the tunnel thermal re- 
sponse and the in-flight response of an electrothermal dei- 
cer. 
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TABLE 1 


Deicer Model Used m Variable Thickness Ice Layer Program 

Sample Runs 


Layer 

Thickness 

Thermal 

Thermal 


(in. ) 

Conductivity 
( BTU/hr-ft- °F) 

Di f fusivi ty 
( ft2/hr) 

Aluminum D-Spar 

0.0870 

66.5 

1.65 

Epoxy/Glass 

0.0200 

0.220 

0.00870 

Insulation 

Point Heater 

0.0000 

— 

— 

Epoxy/Glass 

0.0200 

0.220 

0.00870 

Insulation 

Stainless Steel 

0.0120 

8.70 

0.150 

Abrasion Shield 

Ice 

Variable 

1.32 

0.0469 


2 

Heater Power Density = 25 W/in 2 

Outer Convection Coefficient = 200 BTU/hr-ft °F 
Inner Convection Coefficient = 10 BTU/hr-ft - °F 

Ambient Temperature = -28 °F 

Initial Temperature = -28 °F 
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TABLE 2 


Physical Construction and Thermal Properties of Materials 
Used in the Helicopter Blade and Electrothermal Deicer 


a. Sides of the Blade (Positions 7 through 11) 


Layer 

Nodes 

Thickness 
(in. ) 

Thermal 
Conductivity 
(BTU/hr-ft- °F) 

Thermal 
Diffusivity 
( ft 2 /hr) 

Abrasion Shield 

Stainless Steel 

6 

0.030 

8.7 

0. 15 

Adhesive 

Epoxy 

4 

0.0168 

0.1 

0.0058 

Insulation 

Epoxy/Glass 

4 

0.0138 

0.22 

0.0087 

Adhesive 

Epoxy 

4 

0.0082 

0.1 

0.0058 

Heater Element 
Copper 

4 

0.0065 

60.0 

1.15 

Adhesive 

Epoxy 

4 

0.0082 

0. 1 

0.0058 

Insulation 

Epoxy/Glass 

4 

0.138 

0-22 

0.0087 

Adhesive 

Epoxy 

2 

0.0082 

0.1 

0.0058 

Blade Skin 

Stainless Steel 

4 

0.02 

8.7 

0.15 

Film Adhesive 
FM 1000 

3 

0.01 

0.1 

0.0058 

Doubler 

Aluminum 

10 

0.05 

102 

2.83 

Film Adhesive 
FM1000 

3 

0.01 

0.1 

0.0058 

D-Spar 

Aluminum 

31 

0.175 

102 

2.83 
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b. Stagnation Point (Position 6) 


Layer ] 

Modes 

Thickness 
(m. ) 

Thermal 
Conductivity 
( BTU/hr- f t- °F ) 

Thermal 
Dif fusivity 
( ft 2/hr ) 

Abrasion Shield 

Stainless Steel 

6 

0.030 

8.7 

0.15 

Adhesive 

Epoxy 

4 

0.0168 

0.1 

0.0058 

Insulation 

Epoxy/Glass 

4 

0.0138 

0.22 

0.0087 

Adhesive 

Epoxy 

2 

0.0082 

0.1 

0.0058 

Heater Element 
Copper 

2 

0.0065 

60.0 

1.15 

Adhesive 

Epoxy 

2 

0.0082 

0.1 

0.0058 

Insulation 

Epoxy/Glass 

4 

0.138 

0.22 

0.0087 

Adhesive 

Epoxy 

2 

0.0082 

0.1 

0.0058 

Blade Skin 

Stainless Steel 

4 

0.02 

8.7 

0.15 

Film Adhesive 
FM 1000 

3 

0.01 

0.1 

0.0058 

Noseblock 

Brass 

38 

0.7 

64.15 

1 . 32 

Film Adhesive 
FM1000 

3 

0.01 

0 . 1 

0.0058 

D-Spar 

Aluminum 

21 

0 . 125 

102 

2.83 
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TABLE 3 


Test Matrix Conditions for the Dry Runs 


Escort 

Voo 

Too 

a 

PI 

P2 

2 P3 

Leading 

(MPH) 

(°F) 

(Deg. ) 

( W/in 

2 ) 

57 

50 

64 

0 

0 

0 

0 

60 

50 

33 

0 

0 

0 

0 

61 

50 

26 

0 

0 

0 

0 

62 

200 

19 

2 

0 

0 

0 

64 

200 

20 

2 

6.5 

6.5 

6.5 

65 

200 

20 

2 

14 

14 

14 

66 

200 

20 

2 

24 

24 

24 

67 

200 

20 

6 

0 

0 

0 

68 

200 

20 

6 

14 

14 

14 

69 

120 

17 

0 

0 

0 

0 

70 

120 

17 

0 

2.9 

2.9 

2.9 

71 

120 

17 

0 

16 

16 

16 

73 

50 

5 

0 

0 

0 

0 

74 

64 

-3 

0 

0 

0 

0 

75 

275 

-4 

0 

0 

0 

0 

76 

275 

-4 

0 

2.9 

2.9 

2.9 

77 

275 

-4 

0 

16 

16 

16 

78 

275 

-4 

0 

16 

16 

16 

79 

275 

-4 

6 

0 

0 

0 

80 

275 

-4 

6 

2.9 

2.9 

2.9 

81 

275 

-4 

6 

16 

16 

16 

82 

275 

-4 

6 

16 

16 

16 

86 

100 

-9 

0 

0 

0 

0 

87 

200 

-10 

6 

8 

8 

8 

88 

200 

-11 

6 

16 

16 

16 

89 

200 

-11 

6 

24 

24 

24 

90 

200 

-12 

2 

24 

24 

24 

91 

200 

-11 

4 

24 

24 

24 

92 

120 

-15 

0 

8 

8 

8 

93 

120 

-15 

0 

16 

16 

16 

94 

120 

- 15 

n 

O A 
4* *X 

O A 
Z. 

O A 

z. *± 

95 

120 

-15 

0 

16 

16 

24 

96 

120 

-15 

0 

8 

8 

24 

97 

120 

-15 

0 

0 

0 

24 

98 

120 

-16 

2 

24 

24 

24 

99 

120 

-15 

4 

24 

24 

24 

100 

120 

-15 

6 

8 

8 

8 

101 

120 

-15 

6 

16 

16 

16 

102 

120 

-15 

6 

24 

24 

24 

103 

120 

-15 

6 

24 

24 

24 

105 

120 

-14 

6 

16 

24 

16 

106 

120 

-15 

6 

16 

24 

16 

107 

120 

-15 

6 

8 

24 

8 

108 

120 

-15 

6 

0 

24 

0 

109 

120 

-16 

6 

8 

8 

8 

110 

120 

-15 

6 

8 

8 

8 


t or/ t off 
(Sec . ) 


30/30 1 
30/30 1 
10/30 1 


10/30 1 


30/30 1 

10/30 1 


30/30 1 
30/30 1 
10/30 4 


30/30 1 
30/30 1 
10/30 4 


30/30 1 

30/30 1 

20/30 1 

20/30 1 

20/30 1 

30/30 1 

30/30 1 

20/30 1 

20/30 1 

20/30 1 

20/30 1 

20/30 1 

20/30 1 

20/30 1 

30/30 1 

20/30 1 

30/30 4 

20/30 1 

20/30 2 

20/30 2 

20/30 2 

30/30 2 

180/100 1 
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TABLE 4 


Test Condition Matrix for the Wet Runs 


Escort 

V® 

Too 

a 

PI 

P2 

oP3 

t /t ^ 


LWC 

( 

d 

Reading 

(MPH) 

( °F) 

(Deg. 

) 

(W/m ) 

UI 1 oil 
(Sec. ) 

( g/m 3 ) 

( ym) 

115 

120 

57 

0 

0 

0 

0 


— 

1 

.07 

16 

.8 

116 

120 

56 

0 

24 

24 

24 

20/30 

1 

1 

.07 

16 

. 8 

120 

120 

61 

0 

24 

24 

24 

10/30 

1 

1 

. 07 

16 

. 8 

121 

120 

53 

4 

0 

0 

0 


- 

1 

.07 

16 

.8 

122 

120 

47 

6 

0 

0 

0 


- 

1 

.07 

16 

.8 

123 

120 

47 

6 

0 

0 

0 


- 

1 

.02 

11 

. 7 

124 

120 

49 

4 

0 

0 

0 


- 

1 

.02 

11 

. 7 

125 

120 

46 

2 

0 

0 

0 


- 

1 

. 02 

11 

. 7 

126 

120 

46 

0 

0 

0 

0 


- 

1 

. 02 

11 

. 7 

127 

120 

46 

0 

0 

0 

0 


- 

2 

. 49 

35 

. 4 

128 

120 

50 

2 

0 

0 

0 


- 

2 

. 49 

35 

. 4 

129 

120 

46 

4 

0 

0 

0 


- 

2 

. 49 

35 

. 4 

130 

120 

47 

6 

0 

0 

0 


- 

2 

. 49 

35 

. 4 

131 

120 

48 

6 

0 

0 

0 


- 

2 

.24 

23 

. 6 

132 

120 

45 

4 

0 

0 

0 


- 

2 

.24 

23 

. 6 

133 

120 

45 

2 

0 

0 

0 


- 

2 

.24 

23 

. 6 

134 

120 

47 

0 

0 

0 

0 


- 

2 

.24 

23 

. 6 

135 

200 

46 

0 

0 

0 

0 


- 

1 

.20 

18 

.8 

136 

200 

46 

0 

16 

16 

16 

20/30 

1 

1 

.20 

18 

.8 

137 

200 

47 

2 

0 

0 

0 

- 

1 

.20 

18 

. 8 

138 

200 

47 

2 

16 

16 

16 

20/30 

1 

1 

.20 

18 

. 8 

139 

200 

49 

4 

16 

16 

16 

30/30 

1 

1 

.20 

18 

.8 

140 

200 

50 

4 

0 

0 

0 


- 

1 

.20 

18 

. 8 

141 

200 

50 

6 

0 

0 

0 


- 

1 

.20 

18 

.8 

142 

200 

50 

6 

16 

16 

16 

30/30 

1 

1 

.20 

18 

. 8 

143 

200 

49 

6 

0 

0 

0 


- 

1 

. 10 

14 

.8 

144 

200 

46 

4 

0 

0 

0 


- 

1 

. 10 

14 

.8 

145 

200 

45 

2 

0 

0 

0 


- 

1 

. 10 

14 

. 8 

146 

200 

43 

0 

0 

0 

0 


- 

1 

. 10 

14 

. 8 

147 

200 

44 

0 

0 

0 

0 


- 

0 

.96 

18 

.0 

148 

200 

45 

0 

16 

16 

16 

30/30 

1 

0 

.96 

18 

.0 

149 

200 

44 

2 

0 

0 

0 


- 

0 

.96 

18 

.0 

150 

200 

45 

2 

16 

16 

16 

30/30 

1 

0 

.96 

18 

.0 

151 

200 

44 

4 

0 

0 

0 


- 

0. 

.96 

18. 

0 

152 

200 

45 

4 

16 

16 

16 

30/30 

1 

0 

.96 

18. 

.0 

153 

200 

44 

6 

0 

0 

0 


- 

0 

.96 

18. 

.0 

154 

200 

44 

6 

16 

16 

16 

30/30 

1 

0 

.96 

18 

.0 

155 

280 

49 

6 

0 

0 

0 


- 

0 

.36 

11 . 

.0 

156 

280 

49 

4 

0 

0 

0 


- 

0 

. 36 

11 , 

.0 

157 

280 

48 

2 

0 

0 

0 


- 

0. 

.36 

11. 

0 

158 

280 

48 

0 

0 

0 

0 


- 

0. 

.36 

11 . 

.0 

159 

280 

48 

0 

0 

0 

0 


- 

1 . 

. 17 

20. 

0 

160 

280 

48 

2 

0 

0 

0 


- 

1 . 

. 17 

20. 

0 

161 

280 

47 

4 

0 

0 

0 


- 

1. 

. 17 

20. 

0 

162 

280 

49 

6 

0 

0 

0 


- 

1. 

. 17 

20. 

0 
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TABLE 5 


Test Condition Matrix for the Ice Accretion Documentation 

Tests 


Escort 

Vco 

Too 

a 

LWC 

d 

^ ice. 
(min) 

Reading 

(MPH) 

( ° F ) 

(Deg. ) 

(g/m 3 ) 

( Mm) 

163 

200 

34 

0 

1.2 

23.1 

5 

164 

200 

31 

0 

1.2 

23.1 

5 

— 

200 

31 

2 

1.2 

23 . 1 

5 

— 

200 

31 

4 

1.2 

23.1 

5 

— 

280 

31 

0 

1.0 

16.8 

5 

— 

200 

29 

0 

1.2 

23 . 7 

5 

170 

200 

24 

0 

1.2 

23.7 

5 

171 

100 

24 

0 

2.7 

38.5 

5 

172 

100 

24 

2 

2.7 

38.5 

5 

— 

100 

16 

0 

2.2 

19.2 

5 

174 

100 

16 

2 

2.2 

19.2 

5 

175 

170 

5 

0 

1.7 

19.2 

5 

176 

170 

5 

4 

1.7 

19.2 

5 

177 

100 

2 

2 

2.3 

25.6 

5 

178 

100 

4 

4 

2.3 

25.6 

5 

179 

100 

0 

0 

1.8 

15.0 

5 

180 

100 

1 

2 

1.8 

15.0 

5 

181 

250 

-2 

0 

1.2 

20.0 

5 

182 

250 

-2 

2 

1.2 

20.0 

5 

183 

280 

31 

0 

1.0 

20.0 

5 

185 

250 

28 

0 

1.0 

20.0 

5 

186 

200 

24 

0 

1.0 

20.0 

5 

187 

150 

21 

0 

1.0 

20.3 

5 

188 

100 

19 

0 

1.0 

19.9 

5 

189 

280 

-3 

0 

1.0 

15.1 

5 

190 

280 

-3 

2 

1.0 

15.1 

5 

191 

200 

-10 

0 

1.2 

20.0 

5 

192 

200 

-11 

2 

1.2 

20.0 

5 

— 

120 

-11 

0 

1.1 

15.3 

5 

— 

120 

-11 

2 

1.1 

15.3 

5 
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TABLE 6 


Test Condition Matrix for the Deicing Tests 


a. Simultaneous Heater Zone Cycling Test 


Escort 

Reading 

V 00 
(MPH) 

T CO 

( °F) 

a 

(Deg. 

PI 

) 

P2 

( W/ir) 

P3 

> 2 ) 

t on /t of£ n 
(Sec . ) 

LWC 
( g/m 3 ) 

d 

(yin) 

. \ce 
( mm 

196 

200 

31 

0 

8 

8 

8 

60/60 

1 

1.2 

18.4 

5 

197 

200 

31 

0 

16 

16 

16 

10/30 

2 

1.2 

18.4 

5 

198 

200 

31 

2 

8 

8 

8 

10/30 

1 

1.2 

18 4 

5 

199 

200 

31 

2 

16 

16 

16 

10/60 

2 

1.2 

18.4 

5 

200 

200 

31 

2 

8 

8 

8 

10/60 

2 

1.2 

18.4 

5 

201 

200 

31 

4 

8 

8 

8 

10/60 

1 

1.2 

18.4 

5 

202 

200 

31 

4 

16 

16 

16 

10/60 

1 

1.2 

18.4 

5 

204 

280 

31 

0 

16 

16 

16 

10/60 

1 

1.0 

14 . 5 

5 

205 

280 

31 

0 

8 

8 

8 

10/60 

1 

1.0 

14 . 5 

5 

206 

280 

31 

0 

8 

8 

8 

10/60 

1 

1.0 

14 . 5 

5 

207 

280 

31 

0 

16 

16 

16 

10/60 

1 

1.0 

14 . 5 

5 

209 

250 

28 

0 

16 

16 

16 

10/60 

3 

1.0 

20.0 

5 

210 

250 

28 

0 

8 

8 

8 

10/60 

3 

1.0 

20.0 

5 

211 

200 

24 

0 

8 

8 

8 

10/60 

3 

1.0 

20.0 

5 

212 

200 

24 

0 

8 

8 

8 

20/60 

3 

1.0 

20.0 

5 

213 

200 

24 

0 

16 

16 

16 

20/60 

3 

1.0 

20.0 

5 

214 

150 

21 

0 

16 

16 

16 

20/60 

3 

1.0 

20.3 

5 

215 

150 

21 

0 

16 

16 

16 

20/60 

3 

1.0 

20.3 

5 

216 

100 

20 

0 

16 

16 

16 

20/60 

3 

1.0 

19.9 

5 

217 

100 

20 

0 

24 

24 

24 

10/60 

3 

1.0 

19.9 

5 

218 

170 

5 

0 

24 

24 

24 

10/60 

3 

1.7 

23 . 1 

5 

219 

170 

6 

0 

24 

24 

24 

10/60 

3 

1 . 7 

23 . 1 

2.5 

223 

200 

29 

0 

8 

8 

8 

10/60 

3 

1.2 

19.0 

2.5 

225 

200 

29 

0 

16 

16 

16 

10/60 

3 

1.2 

19.0 

2.5 

226 

200 

24 

0 

8 

8 

8 

10/30 

4 

1.2 

19.0 

2.5 

227 

200 

23 

0 

16 

16 

16 

10/30 

3 

1.2 

19.0 

2.5 

228 

100 

24 

0 

8 

8 

8 

10/30 

3 

2.7 

29.3 

2.5 

229 

100 

24 

0 

16 

16 

16 

10/30 

3 

2.7 

29.3 

2.5 

230 

100 

25 

0 

24 

24 

24 

10/30 

3 

2.7 

29.3 

2.5 

231 

100 

25 

2 

8 

8 

8 

10/30 

3 

2 . 7 

29.3 

2.5 

232 

100 

24 

2 

16 

16 

16 

20/40 

3 

2.7 

29.3 

2.5 

233 

100 

24 

2 

24 

24 

24 

10/40 

3 

2 . 7 

29.3 

2.5 

234 

100 

16 

0 

8 

8 

8 

20,30/60 

3 

2.2 

19.2 

2.5 

238 

200 

29 

0 

16 

16 

16 

10/60 

3 

1.2 

19.0 

5 

239 

200 

29 

0 

8 

8 

16 

10/60 

3 

1.2 

19.0 

5 

240 

200 

29 

0 

0 

0 

16 

10/60 

3 

1.2 

19.0 

5 

241 

100 

16 

0 

16 

16 

16 

20/60 

3 

2.2 

19.2 

2 . 5 

242 

100 

16 

0 

24 

24 

24 

10/60 

3 

2.2 

19.2 

2 . 5 

243 

100 

16 

2 

24 

24 

24 

10/60 

3 

2.2 

19.2 

2.5 

244 

100 

16 

2 

16 

16 

16 

20/60 

4 

2.2 

19.2 

2.5 

245 

100 

16 

2 

8 

8 

8 

30/60 

4 

2 . 2 

19.2 

2.5 

246 

100 

-1 

2 

8 

8 

8 

30/60 

4 

1.8 

15.0 

2.8 

247 

100 

-1 

2 

16 

16 

16 

20/60 

4 

1.8 

15.0 

2.5 

251 

170 

5 

0 

24 

24 

24 

10/60 

4 

1 . 7 

19.7 

2.5 
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252 

170 

5 

0 

16 

16 

16 

20/60 

4 

1.7 

19.7 

2.5 

253 

170 

5 

0 

8 

8 

8 

30/30 

4 

1.7 

19.7 

2.5 

254 

170 

5 

0 

16 

16 

24 

10/60 

4 

1.7 

19.7 

2.5 

255 

170 

5 

0 

8 

8 

24 

10/60 

2 

1.7 

19.7 

2.5 

256 

170 

5 

4 

24 

24 

24 

10/60 

2 

1.7 

19.7 

2.5 

257 

170 

5 

4 

16 

16 

16 

30/30 

4 

1.7 

19.7 

2.5 

258 

170 

5 

4 

8 

8 

8 

60/30 

4 

1.7 

19.7 

2.5 

259 

100 

2 

2 

8 

8 

8 

60/30 

3 

2.3 

25.6 

2.5 

260 

100 

2 

2 

16 

16 

16 

30/30 

4 

2.3 

25.6 

2.5 

261 

100 

2 

2 

24 

24 

24 

10/30 

4 

2.3 

25.6 

2.5 

262 

100 

2 

4 

24 

24 

24 

10/30 

4 

2.3 

25.6 

2.5 

263 

100 

2 

4 

16 

16 

16 

30/30 

4 

2.3 

25.6 

2.5 

264 

100 

2 

4 

8 

8 

8 

60/30 

4 

2.3 

25.6 

2 . 5 

265 

100 

-1 

0 

8 

8 

8 

60/30 

4 

2.3 

25.6 

2.5 

266 

100 

-1 

0 

16 

16 

16 

30/30 

4 

2.3 

25.6 

2.5 

267 

100 

-1 

0 

24 

24 

24 

10/30 

4 

2.3 

25.6 

2.5 

269 

100 

-1 

2 

8 

8 

8 

60/30 

3 

1.8 

15.0 

2.5 

270 

100 

-1 

2 

16 

16 

16 

30/30 

3 

1.8 

15.0 

2.5 

271 

100 

-1 

2 

24 

24 

24 

10/30 

3 

1.8 

15.0 

2.5 

272 

250 

-1 

0 

8 

8 

8 

60/30 

3 

1.2 

20.0 

2.5 

274 

250 

-2 

0 

16 

16 

16 

30/30 

4 

1.2 

20.0 

2.5 

275 

250 

-3 

0 

24 

24 

24 

10/30 

4 

1.2 

20.0 

2.5 

277 

250 

-3 

2 

8 

8 

8 

60/30 

3 

1.2 

20.0 

0.5 

278 

250 

-2 

2 

16 

16 

16 

30/30 

3 

1.2 

20.0 

0.5 

280 

250 

-2 

2 

24 

24 

24 

10/10 

8 

1-2 

20.0 

0.5 

281 

280 

-3 

0 

8 

8 

8 

60/30 

3 

1.0 

15.1 

2.5 

282 

280 

-4 

0 

16 

16 

16 

20/10 

10 

1.0 

15.1 

0.5 

283 

280 

-4 

0 

24 

24 

24 

10/10 

12 

1.0 

15.1 

— 



b 

. Phased Heater 

Cycling 

Test 





Escort 

Reading 

v„ 

(MPH) 

Too 

( °F) 

a 

(Deg 

P 9 HC 
. ) ( W/in j 

ox/ t off 
(Sec . ) 

n 

LWC 

( g/m 

d 

(ym) 

(mifil 

286 

200 

29 

0 

8 

1 

10/60 

3 

1.2 

19.0 

2.5 

287 

200 

29 

0 

8 

1 

10/60 

1 

1.2 

19.0 

2.5 

288 

100 

24 

0 

8 

1 

10/60 

3 

2.7 

29.3 

2.5 

289 

100 

24 

0 

8 

▲ 

10/60 

3 

2 . / 

29 . 3 

2 . 5 

290 

100 

24 

0 

16 

1 

10/60 

3 

2.7 

29.3 

2 . 5 

291 

100 

16 

0 

16 

1 

10/60 

3 

2.2 

19.2 

2 . 5 

292 

100 

16 

0 

8 

1 

10/60 

3 

2.2 

19.2 

2 . 5 

293 

100 

-1 

0 

16 

1 

10/60 

3 

1.8 

15.0 

2.5 

294 

100 

-1 

0 

16 

1 

10/60 

3 

1.8 

15.0 

2.5 

295 

200 

-11 

0 

16 

1 

10/60 

3 

1.2 

20.0 

2.5 

298 

200 

29 

0 

8 

2 

10/40 

3 

1.2 

19.0 

2.5 

300 

200 

29 

0 

8 

2 

10/40 

3 

1.2 

19.0 

5.0 

301 

200 

24 

0 

16 

1 

10/60 

3 

1.0 

20.0 

2.5 

302 

200 

24 

0 

24 

1 

10/60 

3 

1.0 

20.0 

2.5 

303 

200 

24 

0 

8 

2 

10/40 

3 

1.0 

20.0 

2.5 

304 

200 

24 

0 

16 

2 

10/40 

3 

1.0 

20.0 

2.5 

305 

100 

16 

0 

16 

2 

10/40 

3 

2.2 

19.2 

2.5 

306 

100 

16 

0 

24 

3 

20/40 

3 

2.2 

19.2 

2.5 

307 

100 

16 

0 

16 

3 

20/40 

3 

2.2 

19 . 2 

2 . 5 
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308 

100 

16 

0 

24 

3 

20/40 

3 

2.2 

19.2 

2.5 

309 

170 

5 

0 

24 

3 

20/40 

4 

1.7 

19.7 

2.5 

310 

100 

0 

0 

24 

3 

20/40 

4 

1.8 

15.0 

2 . 5 

311 

100 

0 

0 

24 

3 

20/40 

4 

1.8 

15.0 

5 

312 

100 

0 

2 

24 

3 

20/40 

4 

1.8 

15.0 

5 

313 

200 

-11 

2 

24 

3 

20/40 

3 

1.2 

20.0 

2.5 

314 

200 

-13 

0 

24 

3 

20/40 

3 

1.2 

20.0 

2.5 

315 

200 

-11 

0 

24 

4 


3 

1.2 

20.0 

2.5 

P = 

Nominal 

Heater 

Power 

Density 








HC = Heater Configuration 

There were four heater configurations. See Figure 50. 
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TABLE 7 


Test Conditions for those Readings Chosen for the Numerical 
and Experimental Data Comparisons 


Escort 

Reading 

Voo 

(MPH) 

o 

8 

a 

(Deg. ) 

P 2 
(W/m z 

t /t 
. on o 

' (Sec 

ff" 

.) 

LWC 
( g/m 

d 

(Mm) 

. ice 
(min) 

70 

120 

17 

0 

2.9 

30/30 

1 

— 

— 

— 

76 

275 

-4 

0 

2.9 

30/30 

1 

— 

— 

— 

92 

120 

-15 

0 

8 

30/30 

1 

— 

— 

— 

93 

120 

-15 

0 

16 

30/30 

1 

— 

— 

— 

94 

120 

-15 

0 

24 

20/30 

1 

— 

— 

— 

197 

200 

31 

0 

16 

10/30 

2 

1.2 

18.4 

5 

209 

250 

28 

0 

16 

10/60 

3 

1.0 

20.0 

5 

213 

200 

24 

0 

16 

20/60 

3 

1.0 

20.0 

5 

234 

100 

16 

0 

8 20 

,30/60 

3 

2.2 

19.2 

2 . 5 

275 

250 

-3 

0 

24 

10/30 

3 

1.2 

20.0 

2.5 


TABLE 8 

Presentation of the Thermocouples Chosen as the 
Representative Set 


Position 

Thermocouple 

Around 

Number 

Blade 

AB H S 


1 

49 



2 

50 

— 

— 

3 

62 

— 

— 

4 

63 

25 

9 

5 

53 

26 

10 

6 

65 

27 

11 

7 

55 

28 

T ^ 
14 

8 

56 

X 

13 

9 

57 

X 

14 

10 

58 

39 

15 

11 

X 

32 

16 


- Indicates that a thermocouple is not present at the 
specified location. 

x Indicates that all of the thermocouples at the specified 
location were mnacurate or discarded. 

AB = Abrasion Shield 
H = Heater 
S = D-Spar 
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TABLE 9 


Presentation of the Zonal Dependent Heater Power Densities 


a. Dry Runs 


Position 

Heater 


Power 

Density 

(W/m 2 ) 


Around 

zone 


Reading Number 


Blade 

1 


70 

76 

92 

93 

94 

2 

3 

4 

8 

2.87 

2 . 87 

7.90 

15 . 8 

23.8 

5 

7 

2 . 87 

2 . 87 

7.91 

15.9 

23.8 

6 

6 

2.86 

2 . 86 

7.87 

15.8 

23 . 7 

7 

5 

2.89 

2.89 

7.97 

16.0 

24.0 

8 

4 

2.96 

2.96 

8.14 

16.3 

24.5 

9 

3 

2.91 

2.91 

8.03 

16.1 

24.2 

10 

2 

2.97 

2.97 

8.20 

16.4 

24.7 

11 

1 

2.97 

2.97 

8.20 

16.4 

24.7 




b . 

Deicing 

Runs 



Position 

Heater 


Power 

Density 

(W/m 2 ) 


Around 

zone 


Reading Number 


Blade 

1 


197 

209 

213 

234 

275 

2 

3 

4 

8 

15.8 

15.8 

15.8 

7.90 

23.8 

5 

7 

15.9 

15.9 

15.9 

7.91 

23.8 

6 

6 

15.8 

15.8 

15.8 

7.87 

23.7 

7 

5 

16.0 

16.0 

16.0 

7.97 

24.0 

8 

4 

16.3 

16.3 

16.3 

8.14 

24.5 

9 

3 

16.1 

16.1 

16.1 

8.03 

24.2 

10 

2 

16.4 

16.4 

16.4 

8.20 

24.7 

11 

1 

16.4 

16.4 

16.4 

8.20 

24.7 
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Figure 2 . Method of Approximating the Variable 
Thickness Ice Layer 


141 


1 ,k+l 


Figure 3. 



a. Interior Mesh Point at the 
Center of the Gap 



b. Interior Mesh Point at the 
Center of the Heater 


Finite Difference Grid at Selected Points 
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Layer 1 


Figure 4 



Interface 

Ambient 


a. Substrate-Ambient Interface 



Ambient 
Interface 
Abrasion Shield 


b. Abrasion Shield-Ambient 
Interface 


. Finite Difference Grid at Selected Points 
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a. 


Type 1 Element 


b. Type 2 Element 




31 

ay 



It -T ~ 

V 1 ‘ " 



e. Type 5 Element 



f. Type 6 Element 


Figure 5, Nodal Boundary Configurations for the Variable Thickness Ice Layer 
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+ + + + 

® ® ® Ice 

Interface 
Abrasion Shield 

c. Type 9 Element 

Figure 6. Nodal Boundary Configurations for the Variable 
Thickness Ice Layer 
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T i ,k+l ,6 



h - " H 


a. Interior Ice Node 



T i,k-1, 5 


|^— ■■ | 

b. Ice-Abrasion Shield Interface 


Figure ?• The Nodal Energy Balance 
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ration 

















Figure 9. Flow Chart for the Enthalpy Method 
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Deicing Time Predictions Using Marano 




Positive X-Direction (Inches) 


Figure II. Progression of Melt Interface for Various Ice Shapes 
With a Mean Ice Thickness of 0.125 (1/8) Inches 
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0.00 0.025 0.050 0.075 0.100 0.125 

Positive X-Direction (Inches) 

Figure 12. Progression of Melt Interface for Various Ice 

Shapes with a Mean Ice Thickness of 0.0625 (1/16) Inches 
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Notei All Dimensions in Inches 
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Details of a UH-1H Helicopter Blade Construction 
(All Dimensions in Inches) 



Abrasion Shield 


Insulation 


Insulation 


Blade Skin 



Epoxy Glue 
Copper Heater 
Epoxy Glue-— 


Abrasion Shield 




Blade Skin 



FM 100 

Film Adhesive 


Doubler 



D-Spar 
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Figure 16. Thermocouple location and Numbering for Thermocouples 
on the Inner Side of the D-Spar 





Top of Blade 



Figure 18. Thermocouple location and Numbering for Thermocouples 
at the Outer Side of the Abrasion Shield 








for Escort Reading 70 


1 1 3 


o 



«) O O O O O 

O • 4 # 

u> 3&ni va3dwai 
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Figure 19. (Reading 70 Continued) 



Abrasion Shield 1VJ L O Abrasion Shield 
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Figure 19* (Reading ?0 Continued 



(i > jS(U*83d*31 
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19. (Reading 70 





Abrasion Shield 
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Figure 19. (Reading 70 Continued) 


Abrasion Shield 
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figure 20. Response of Selected Thermocouples 
for Escort Reading ?6 



Abrasion Shield ._ u O Abrasion Shield 
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Figure 20 . (Reading 76 Continued) 






Q Abrasion Shield ,J r Q Abrasion Shield 

i 

C] Heater ,J CH Heater 




u) 3aniva3dM3i 


On 

C 

o 


<0 

(2 



168 


Figure 20. (Reading ? 6 Continued) 



O Abrasion Shield •« - O Heater 



U> 3tirnv«3dHJl 
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Figure 20. (Reading ?6 Continued) 



Figure 21. Response of Selected Thermocouples 
for Escort Reading 92 


Q Abrasion Shield >**h Q Abrasion Shield 
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Figure 21. (Reading 92 Continued) 


Q Abrasion Shield Q Abrasion Shield 
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Figure 21. (Reading 92 Continued) 




< J> 3Bn)VB3dH3l 
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Figure 21. (Reading 92 Continued) 







Figure 22. Response of Selected Thermocouples 
for Escort Reading 93 



Abrasion Shield >•<>[- Q Abrasion Shield 
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Figure 22. (Reading 93 Continued) 
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Figure 22. (Reading 93 3 




O Abrasion Shield - O Abrasion Shield 



(j> 3aniv«3dM3i 
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Figure 22. (Reading 93 Continued) 



1*0 



<J> J«ru»a3dM3i 
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Figure 22. (Reading 93 



Figure 23. Response of Selected Thermocouples 
for Escort Reading 9^ 






Abrasion Shield >«(- , O Abrasion Shield 
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Figure 23. (Reading 9*+ 


Q Abrasion Shield - Q Abrasion Shield 
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Figure 23. (Reading 9^ Continued) 






Q Abrasion Shield 



Figure 23. (Reading 9 ^ Continued) 



Abrasion Shield 


i 
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Figure 24. Response of Selected Thermocouples 
for Escort Reading 197 



Abrasion Shield " 3 [ O Abrasion shield 
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Figure 24. (Reading 197 Continued) 





Q Abrasion Shield “ 0 - O Abrasion Shield 

I I Heater >•" - CD Heater 
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Figure 24. (Reading 197 Continued) 


O Abrasion Shield <« - Q Abrasion Shield 
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Figure 24. (Reading 197 Continued) 
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Figure 24. (Reading 197 3on 
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Figure 25. Response of Selected Thermocouples 
for Escort Reading 209 


o Abrasion Shield o Abrasion Shield 
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Figure 25. (Reading 209 Continued) 
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Figure 25. (Reading 209 Continued) 
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Figure 25. (Reading 209 Continued) 









Figure 26. Response of Selected Thermocouples 
for Escort Reading 213 


o Abrasion Shield “ c j" O Abrasion Shield 


t j > janivaiaHJi 
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Figure 26. (Reading 213 Continued) 




Figure 26. (Reading 213 Continued) 


Abrasion Shield ! O Abrasion Shield 



Figure 26. (Reading 213 Continued) 


Abrasion Shield 



Figure 26. (Reading 213 Continued) 






Figure 26. (Reading 213 Continued) 




Abrasion Shield 



Figure 27. Response of Selected Thermocouples 
for Escort Reading 23^- 


Abrasion Shield 1,0 r Q Abrasion Shield 
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Figure 2?. (Reading 234 Continued) 
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Figure 27* (Reading 23^ G' 
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Figure 27. (Reading 234 Continued) 



Q Abrasion Shield t«L- Q Abrasion Shield 
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Figure 2?. (Reading 23^ Continued) 



Abrasion Shield 
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Figure 2?. (Reading 23^+ 
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Figure 28. Response of Selected Thermocouples 
for Escort Reading 2?5 


Abrasion Shield Q Abrasion Shield 
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Figure 28 . (Reading 275 Continued) 
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Figure 28. (Reading 275 Continued) 




O Abrasion Shield O Abrasion Shield 
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Figure 28* (Reading 275 Continued) 
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Figure 28 . (Reading 275 Continued) 
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Figure 31. Presentation of Positional Trends at the 
Abrasion Shield for Reading 19? 
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Figure 32. Presentation of Positional Trends at the Abrasion 
Shield for Reading 209 
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First 

Peak 


Second 








Position 


Figure 35. Presentation of Positional Trends at the 
Abrasion Shield for Reading 2 75 








12 34 5 6 7 8 9 10 11 

Position 


O Average of all three cycles and all five deicing runs. 


Figure 36. Presentation of Positional Trends at the 
Heater for the Deicing Runs 
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Figure 38- Numerical and Experimental Comparison for Reading 70 
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Figure J8. (Reading 70 Ot 
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Figure 39. Numerical and Experimental Comparison 
for Reading 76 
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Figure 39. (Reading 76 Continued) 
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Figure 40. (Reading 92 Continued) 
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Figure 41. Numerical and Experimental Comparison 
for Reading 93 
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Figure 41. (Reading 93 Continued) 
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Figure 41. (Reading 93 Continued) 
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Figure U2. Numerical and Experimental Comparison 
for Reading 9^ 
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Figure 42. (Reading 94 Continued) 
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Figure 42. (Reading 94 Continued) 
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Figure 43. Presentation of the Capton Effect 
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Figure 45. Numerical and Experimental r Jomparison 
for Reading 209 
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Figure 46. Numerical and Experimental Comparison 
for Reading 213 
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Figure 47. Numerical and Experimental Comparison 
for Reading 23^ 
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Figure 4?. (Reading 234 Continued) 



Figure 4?. (Reading 23 4 Continued) 
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Figure 47« (Reading 23^- Continued) 
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Figure 48. Numerical and Experimental Comparison 
for Reading 275 
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Figure 49* Response of Thermocouples at Positions 
6 and 7 for Reading 234 
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Appendix A 


The Enthalpy Method simulation of a phase change process 
has been shown to give reliable results, as documented m 
Marano [ 1A ] . However, due to the finite distance between 
nodes m a finite-difference formulation of a problem, the 
numerical results for the temperature versus time behavior 
of a node that undergoes a phase change may show an unreal- 
istic response after melting begins. This will occur unless 
a very large number of nodes are used in the ice layer. 

As illustrated in Figure lAa (Figure 12 of ( 1A ] ) for the 
abrasion shield ice interface temperature response for a 
deicer pad, the temperature, after it reaches 32° F, remains 
at 32 °F for a certain time and then oscillates after com- 
plete melting occurs. Both the plateau and the oscillation 
frequency are nodal dependent . This can be seen by compar- 
ing the two curves for 20 and 60 nodes m the ice layer. 
These oscillations are attributed to the fact that a node in 
the ice layer remains at the melting point for a finite 
period of time until the volume of ice represented by the 
node has entirely melted. 

Voller and Cross (2A,3A) have shown, by comparing analyt- 
ical and numerical solutions to simple phase change problems 


257 


using the Enthalpy Method, that the numerical solutions os- 
cillate around the true solutions. They have also developed 
a criterion for determining the points of correspondence be- 
tween the true and oscillating solutions. By finding these 
points of correspondence, accurate response curves can be 
obtained. 

In Essence, Voller and Cross have shown that the liquid- 
solid interface is exactly located at the center of a node 
when the nodal enthalpy is Hsmp + (Hlmp - Hsmp)/2. For the 
ice-water system, this value is [925.63 + (9032.31)/2] = 
5441.79 BTU/ft"^ . By plotting the abrasion shield surface 
temperature at these specific times, the true solution is 
obtained. This procedure was used to replot the 20 and 60 
node curves in Figure lAa, as well as to plot data for 30, 
40 and 90 nodes in the ice layer. The result is shown m 
Figure lAb (Figure 13 of [ 1A ] ) 

To be more specific, the 90 node curve in Figure lAb was 
obtained in the following manner. At 4.9 sec., the inter- 
face node begins to melt. At 7.23 sec., the enthalpy for 

3 

the second node reaches 5441.79 BTU/ft , meaning that the 
ice-water interface is exactly at the second node. The 
first triangle m Figure lAb was taken at this time, when 
the abrasion shield surface temperature was 36.66 °F. At 
9.13 sec., the enthalpy of the third node reaches this same 
value, and the surface temperature at that time is plotted 
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as the second triangle. This procedure is repeated as time 
increases. As stated in Marano [lAj, 30 nodes per 0.25 
inches of ice was found to be the practical minimum number 
of nodes for sufficient accuracy. 

Finally, as stated above, this temperature plateau and 
the succeeding oscillations can be completely elimmatecj by 
using a very large number of nodes m the ice layer. This 
would enable the direct program output to be plotted instead 
of following the procedure described m the preceedmcj para- 


graph . 
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